TABLE OF CONTENTS
- ABINIT/berryphase_new
- ABINIT/init_e_field_vars
- ABINIT/initberry
- ABINIT/m_berryphase_new
- ABINIT/prtefield
- ABINIT/update_e_field_vars
ABINIT/berryphase_new [ 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.
SOURCE
159 subroutine berryphase_new(atindx1,cg,cg13,cprj,dtefield,dtfil,dtset,psps,& 160 & gprimd,hdr,indlmn,kg,lmnmax,mband,mcg,mcg13,mcprj,& 161 & mkmem,mpi_enreg,mpw,my_natom,natom,npwarr,nsppol,ntypat,& 162 & nkpt,calc_pol_ddk,pawrhoij,pawtab,pel,pelev,pion,ptot,red_ptot,pwind,& !!REC 163 & pwind_alloc,pwnsfac,& 164 & rprimd,save_cg13,typat,ucvol,unit_out,usecprj,usepaw,xred,zion) 165 166 !Arguments ------------------------------------ 167 integer, intent(in) :: lmnmax,mband,mcg,mcg13,mcprj,mkmem,mpw,my_natom,natom,nkpt 168 integer, intent(in) :: nsppol,ntypat,calc_pol_ddk 169 integer, intent(in) :: pwind_alloc,unit_out,usecprj,usepaw 170 real(dp), intent(in) :: ucvol 171 logical,intent(in) :: save_cg13 172 type(MPI_type), intent(in) :: mpi_enreg 173 type(datafiles_type), intent(in) :: dtfil 174 type(dataset_type), intent(in) :: dtset 175 type(pseudopotential_type),intent(in) :: psps 176 type(efield_type), intent(inout) :: dtefield 177 type(hdr_type), intent(inout) :: hdr 178 !arrays 179 integer, intent(in) :: atindx1(natom),indlmn(6,lmnmax,ntypat),kg(3,mpw*mkmem) 180 integer, intent(in) :: npwarr(nkpt),pwind(pwind_alloc,2,3) 181 integer, intent(in) :: typat(natom) 182 real(dp), intent(in) :: cg(2,mcg),gprimd(3,3) 183 real(dp), intent(in) :: pwnsfac(2,pwind_alloc) 184 real(dp), intent(in) :: rprimd(3,3),zion(ntypat) 185 real(dp), intent(inout) :: xred(3,natom) 186 real(dp), intent(out) :: pel(3),pelev(3),pion(3) 187 real(dp), intent(out) :: ptot(3),red_ptot(3) !!REC 188 real(dp),intent(out) :: cg13(2,mcg13,3) 189 type(pawrhoij_type), intent(in) :: pawrhoij(my_natom*usepaw) 190 type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw) 191 type(pawcprj_type),intent(in) :: cprj(natom,mcprj*usecprj) 192 193 !Local variables ------------------------- 194 integer :: count,count1,dest,fdir,unt 195 integer :: iatom,iband,icg,icg1,idir,idum,ikpt1i_sp 196 integer :: ierr,ifor,ikg,ikpt,ikpt1,ikpt_loc!,ikpt2,ikpt2i,npw_k2, itrs 197 integer :: icp1, icp2,icpgr_offset,iproc 198 ! integer :: ii ! appears commented out below in a debug section 199 integer :: inibz,ikpt1i 200 integer :: isppol,istr,itypat,jband,jkpt,jkstr,jsppol 201 integer :: det_inv_smat, det_smat, inv_smat 202 integer :: maxbd,mcg1_k 203 integer :: minbd,my_nspinor,nband_k,ncpgr,nfor,npw_k1,ntotcp,n2dim,nproc,pertcase 204 integer :: response,shiftbd,source,spaceComm,tag 205 integer :: jj,jstr,kk,ineigh_str 206 integer :: istep,jstep,kpt_mark(dtefield%fnkpt),nkstr,nstr,iunmark,berrystep 207 integer :: jkpt2, jkpt2i 208 real(dp) :: det_mod,dkinv,dphase,dtm_real,dtm_imag,fac,gmod,phase0 209 real(dp) :: pol,polbtot,polion,politot,poltot,rho 210 logical :: calc_epaw3_force,calc_epaw3_stress,efield_flag 211 integer :: polflag, ddkflag 212 213 !!REC start 214 integer :: jump 215 real(dp),save ::pol0(3) 216 logical, save :: first=.true. 217 logical :: lexist 218 !!REC end 219 real(dp) :: dphase_new,dphase_init 220 character(len=fnlen) :: fiwf1o 221 character(len=500) :: message 222 type(wvl_wf_type) :: wfs 223 type(wvl_internal_type) :: wvl 224 integer,allocatable :: dimlmn(:),ikpt1_recv(:), sflag_k(:)!,pwind_k(:) 225 integer,allocatable :: ikpt3(:), ikpt3i(:), sflag_k_mult(:,:),nattyp_dum(:),npw_k3(:) 226 integer,allocatable :: idxkstr_mult(:,:), pwind_k_mult(:,:),itrs_mult(:) 227 real(dp) :: det_average(2),dk(3),dtm_k(2),gpard(3),pel_cart(3),pion_cart(3) 228 real(dp) :: polb(nsppol),ptot_cart(3),rel_string(2),xcart(3,natom) 229 real(dp) :: delta_str(2),dist_,dstr(2) 230 real(dp),allocatable :: buffer(:,:),buffer1(:),buffer2(:) 231 real(dp),allocatable :: cg1(:,:),cg1_k(:,:),cgq(:,:) 232 real(dp),allocatable :: det_string(:,:),dsdr(:,:,:,:,:),dsdr_sum(:,:,:),dsds_sum(:,:,:) 233 real(dp),allocatable :: eig_dum(:),epawf3_str(:,:,:),epaws3_str(:,:,:) 234 real(dp),allocatable :: occ_dum(:),polberry(:),resid(:),pwnsfac_k(:,:) 235 real(dp),allocatable :: smat_inv(:,:,:),smat_k(:,:,:),smat_k_paw(:,:,:) 236 real(dp),allocatable :: str_flag(:) 237 ! real(dp),allocatable :: dist_str(:,:),det_string_test(:,:) 238 real(dp),allocatable :: dtm_mult(:,:,:), coef(:,:), polb_mult(:,:) 239 type(pawcprj_type),allocatable :: cprj_k(:,:),cprj_kb(:,:),cprj_buf(:,:),cprj_gat(:,:) 240 type(pawcprj_type),allocatable :: cprj_fkn(:,:),cprj_ikn(:,:) 241 242 ! integer :: bband,bbs,bra_start,bra_end,ilmn,ipw,ispinor,jlmn,kband,kbs,ket_start,ket_end,klmn,npw_k 243 ! integer :: nspinor,spnipw,spnshft 244 ! real(dp) :: err_ovlp,mag_ovlp,max_err_ovlp, ovlp_r, ovlp_i, paw_r, paw_i 245 ! real(dp) :: tot_r, tot_i 246 ! real(dp),allocatable :: bra(:,:),ket(:,:) 247 ! complex(dpc) :: cpb,cpk,cterm 248 249 !BEGIN TF_CHANGES 250 integer :: me 251 !END TF_CHANGES 252 253 !no_abirules 254 255 ! *********************************************************************** 256 257 !DEBUG 258 !write(std_out,*)' berryphase_new : enter' 259 !do ii=1,pwind_alloc 260 !write(std_out,*)ii,pwnsfac(:,ii) 261 !end do 262 !stop 263 !ENDDEBUG 264 265 !Init MPI 266 spaceComm=mpi_enreg%comm_cell 267 nproc=xmpi_comm_size(spaceComm) 268 me=mpi_enreg%me_kpt 269 270 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 271 272 polflag = 1 273 ddkflag = 1 274 if (calc_pol_ddk == 1) then 275 ddkflag = 0 276 else if (calc_pol_ddk == 2) then 277 polflag = 0 278 end if 279 280 !allocate(pwind_k(mpw)) 281 ABI_MALLOC(pwnsfac_k,(4,mpw)) 282 ABI_MALLOC(sflag_k,(dtefield%mband_occ)) 283 !pwind_k(:) = 0 284 pwnsfac_k(1,:) = 1.0_dp ! bra real 285 pwnsfac_k(2,:) = 0.0_dp ! bra imag 286 pwnsfac_k(3,:) = 1.0_dp ! ket real 287 pwnsfac_k(4,:) = 0.0_dp ! ket imag 288 289 if (maxval(dtset%istwfk(:)) /= 1) then 290 write(message, '(a,a,a)' )& 291 & 'This routine does not work yet with istwfk /= 1.',ch10,& 292 & 'This should have been tested previously ...' 293 ABI_BUG(message) 294 end if 295 296 if (usepaw == 1 .and. usecprj /= 1) then 297 message = ' PAW calculation but cprj datastructure has not been allocated !' 298 ABI_BUG(message) 299 end if 300 301 if (save_cg13 .AND. (ddkflag /= 1) ) then 302 message = ' cg13 output requested but ddkflag not set ' 303 ABI_BUG(message) 304 end if 305 306 if (save_cg13 .AND. (mcg13 /= mcg) ) then 307 message = ' cg13 output requested but mcg13 /= mcg ' 308 ABI_BUG(message) 309 end if 310 311 if (save_cg13 .AND. (size(cg13) /= 2*mcg13*3 ) ) then 312 message = ' cg13 output requested but cg13 size incorrect ' 313 ABI_BUG(message) 314 end if 315 316 ! useful flags for various efield possibilities 317 efield_flag = ( dtset%berryopt == 4 .or. dtset%berryopt == 6 .or. dtset%berryopt == 7 .or. & 318 & dtset%berryopt ==14 .or. dtset%berryopt ==16 .or. dtset%berryopt ==17 ) 319 calc_epaw3_force = ( efield_flag .and. dtset%optforces /= 0 .and. usepaw == 1 ) 320 calc_epaw3_stress = ( efield_flag .and. dtset%optstress /= 0 .and. usepaw == 1 ) 321 322 mcg1_k = mpw*mband 323 shiftbd = 1 324 if (ddkflag==1) then 325 ABI_MALLOC(cg1,(2,mcg)) 326 ABI_MALLOC(eig_dum,(2*mband*mband*nkpt*nsppol)) 327 ABI_MALLOC(occ_dum,(mband*nkpt*nsppol)) 328 eig_dum(:) = zero 329 occ_dum(:) = dtefield%sdeg 330 end if 331 332 !initialize variable tied to multiple step computation 333 berrystep=dtset%berrystep 334 ABI_MALLOC(ikpt3,(berrystep)) 335 ABI_MALLOC(ikpt3i,(berrystep)) 336 ABI_MALLOC(sflag_k_mult,(dtefield%mband_occ,berrystep)) 337 ABI_MALLOC(npw_k3,(berrystep)) 338 ABI_MALLOC(pwind_k_mult,(mpw,berrystep)) 339 ABI_MALLOC(itrs_mult,(berrystep)) 340 ABI_MALLOC(coef,(berrystep,berrystep)) 341 ABI_MALLOC(polb_mult,(nsppol,berrystep)) 342 !coefficient for berryphase computation 343 coef(:,:) = 0.0_dp 344 do jstep = 1, berrystep 345 coef(jstep,1) = 1.d0/real(jstep*jstep,dp) 346 if(jstep/=1)coef(jstep,1)=coef(jstep,1)/real(1-jstep*jstep,dp) 347 end do 348 do istep = 2, berrystep 349 do jstep = 1, berrystep 350 coef(jstep, istep) = real(istep*istep,dp)*coef(jstep,istep-1) 351 if(jstep /= istep)coef(jstep, istep)=coef(jstep,istep)/real(istep*istep-jstep*jstep,dp) 352 end do 353 end do 354 !the berryphase using the strings of steps dk, 2*dk, ..., istep*dk is : 355 !coef(1,istep)*berryphase(dk) + coef(2,istep)*berryphase(2*dk) + ... + coef(istep,istep)*berryphase(istep*dk) 356 !DEBUG 357 !write(std_out,*)'coef, sum coef' 358 !do istep=1,step 359 !write(std_out,*)coef(:,istep), sum(coef(1:istep,istep)) 360 !end do 361 !ENDDEBUG 362 363 !allocate(dtm(2,dtefield%fnkpt*nsppol)) 364 ABI_MALLOC(dtm_mult,(2,dtefield%fnkpt*nsppol,berrystep)) 365 ABI_MALLOC(cg1_k,(2,mcg1_k)) 366 367 if (usepaw == 1) then ! cprj allocation 368 ncpgr = cprj(1,1)%ncpgr 369 if ( calc_epaw3_force ) then 370 ABI_MALLOC(dsdr_sum,(natom,3,dtefield%fnkpt*nsppol)) 371 ABI_MALLOC(epawf3_str,(natom,3,3)) 372 end if 373 if ( calc_epaw3_stress ) then 374 ABI_MALLOC(dsds_sum,(natom,6,dtefield%fnkpt*nsppol)) 375 ABI_MALLOC(epaws3_str,(natom,3,6)) 376 end if 377 ABI_MALLOC(dimlmn,(natom)) 378 call pawcprj_getdim(dimlmn,natom,nattyp_dum,ntypat,typat,pawtab,'R') 379 ABI_MALLOC(cprj_k,(natom,dtefield%nspinor*mband)) 380 ABI_MALLOC(cprj_kb,(natom,dtefield%nspinor*mband)) 381 ABI_MALLOC(cprj_gat,(natom,nproc*dtefield%nspinor*mband)) 382 call pawcprj_alloc(cprj_k,ncpgr,dimlmn) 383 call pawcprj_alloc(cprj_kb,ncpgr,dimlmn) 384 call pawcprj_alloc(cprj_gat,ncpgr,dimlmn) 385 if (dtset%kptopt /= 3) then 386 ABI_MALLOC(cprj_ikn,(natom,dtefield%nspinor*mband)) 387 ABI_MALLOC(cprj_fkn,(natom,dtefield%nspinor*mband)) 388 call pawcprj_alloc(cprj_ikn,ncpgr,dimlmn) 389 call pawcprj_alloc(cprj_fkn,ncpgr,dimlmn) 390 end if 391 392 n2dim = dtefield%nspinor*mband 393 ntotcp = n2dim*SUM(dimlmn(:)) 394 if (nproc>1) then 395 ABI_MALLOC(cprj_buf,(natom,dtefield%nspinor*mband)) 396 call pawcprj_alloc(cprj_buf,ncpgr,dimlmn) 397 end if 398 399 if ( efield_flag ) then 400 write(message,'(2a,i5,2a)')ch10,& 401 & ' nkpt = ',nkpt,ch10,' copy cprj to dtefield%cprj ' 402 call wrtout(std_out,message,'COLL') 403 404 do isppol = 1, nsppol 405 406 ikpt_loc = 0 407 ikpt1 = 0 408 do while (ikpt_loc < mkmem) 409 410 if (ikpt_loc < mkmem) ikpt1 = ikpt1 + 1 411 if ((ikpt1 > nkpt).and.(ikpt_loc < mkmem)) exit 412 nband_k = dtset%nband(ikpt1) 413 414 if ( (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt1,1,nband_k,isppol,me)) .and. & 415 & (ikpt_loc <= mkmem) ) cycle 416 417 ikpt_loc = ikpt_loc + 1 418 419 ABI_MALLOC(ikpt1_recv,(nproc)) 420 call xmpi_allgather(ikpt1,ikpt1_recv,spaceComm,ierr) 421 call pawcprj_get(atindx1,cprj_k,cprj,natom,1,(ikpt_loc-1)*nband_k*my_nspinor,ikpt1,0,isppol,mband,& 422 & mkmem,natom,nband_k,nband_k,my_nspinor,nsppol,0,& 423 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 424 call pawcprj_mpi_allgather(cprj_k,cprj_gat,natom,n2dim,1,dimlmn,ncpgr,nproc,spaceComm,ierr,rank_ordered=.true.) 425 do iproc = 1, nproc 426 icp2=nband_k*(iproc-1)*my_nspinor 427 call pawcprj_get(atindx1,cprj_k,cprj_gat,natom,1,icp2,ikpt1,0,isppol,mband,& 428 & nproc,natom,nband_k,nband_k,my_nspinor,1,0,& 429 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 430 icp1 = nband_k*(ikpt1_recv(iproc)-1)*my_nspinor 431 call pawcprj_put(atindx1,cprj_k,dtefield%cprj,natom,1,icp1,ikpt1,0,isppol,& 432 & mband,dtefield%fnkpt,natom,nband_k,nband_k,dimlmn,my_nspinor,nsppol,0,& 433 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 434 end do 435 ABI_FREE(ikpt1_recv) 436 437 end do ! close loop over k-points 438 439 end do ! end loop over nsppol 440 441 end if ! end check on efield 442 end if 443 444 !!======================================= 445 !! code to test orthonormality of cg_k 446 !!======================================= 447 ! 448 !ikpt = 3 449 !npw_k = npwarr(ikpt) 450 !isppol = 1 451 !nband_k = dtefield%mband_occ 452 !ABI_MALLOC(bra,(2,npw_k*my_nspinor)) 453 !ABI_MALLOC(ket,(2,npw_k*my_nspinor)) 454 !max_err_ovlp=0.0 455 !call pawcprj_get(atindx1,cprj_k,cprj,natom,1,dtefield%cprjindex(ikpt,isppol),ikpt,0,isppol,mband,& 456 !& mkmem,natom,nband_k,nband_k,my_nspinor,nsppol,0) 457 !do bband = 1, nband_k 458 !bra_start = dtefield%cgindex(ikpt,nsppol)+1+(bband-1)*npw_k*my_nspinor 459 !bra_end = bra_start + npw_k*my_nspinor - 1 460 !bra(1:2,1:npw_k*my_nspinor) = cg(1:2,bra_start:bra_end) 461 !do kband = 1, nband_k 462 !ket_start = dtefield%cgindex(ikpt,nsppol)+1+(kband-1)*npw_k*my_nspinor 463 !ket_end = ket_start + npw_k*my_nspinor - 1 464 !ket(1:2,1:npw_k*my_nspinor) = cg(1:2,ket_start:ket_end) 465 ! 466 !tot_r = 0.0; tot_i = 0.0 467 !do ispinor = 1, my_nspinor 468 !ovlp_r = 0.0; ovlp_i = 0.0 469 !spnshft = (ispinor-1)*npw_k 470 !do ipw = 1, npw_k 471 !spnipw = ipw + spnshft 472 !ovlp_r = ovlp_r + bra(1,spnipw)*ket(1,spnipw)+bra(2,spnipw)*ket(2,spnipw) 473 !ovlp_i = ovlp_i - bra(2,spnipw)*ket(1,spnipw)+bra(1,spnipw)*ket(2,spnipw) 474 !end do ! end loop over ipw 475 !paw_r = 0.0; paw_i = 0.0 476 !do iatom = 1, natom 477 !itypat = typat(iatom) 478 !do ilmn = 1, dtefield%lmn_size(itypat) 479 !do jlmn = 1, dtefield%lmn_size(itypat) 480 !klmn=max(ilmn,jlmn)*(max(ilmn,jlmn)-1)/2 + min(ilmn,jlmn) 481 !bbs = my_nspinor*(bband-1)+ispinor 482 !kbs = my_nspinor*(kband-1)+ispinor 483 !cpb=cmplx(cprj_k(iatom,bbs)%cp(1,ilmn),cprj_k(iatom,bbs)%cp(2,ilmn)) 484 !cpk=cmplx(cprj_k(iatom,kbs)%cp(1,jlmn),cprj_k(iatom,kbs)%cp(2,jlmn)) 485 !cterm = conjg(cpb)*pawtab(itypat)%sij(klmn)*cpk 486 !paw_r = paw_r + real(cterm) 487 !paw_i = paw_i + aimag(cterm) 488 !end do ! end loop over jlmn 489 !end do ! end loop over ilmn 490 !end do ! end loop over iatom 491 !tot_r = tot_r + ovlp_r + paw_r 492 !tot_i = tot_i + ovlp_i + paw_i 493 !end do ! end loop over ispinor 494 ! 495 !! write(std_out,'(a,2i4,2es16.8)')' JWZ Debug: berryphase_new bband kband ovlp : ',& 496 !!& bband,kband,tot_r,tot_i 497 !mag_ovlp = tot_r*tot_r + tot_i*tot_i 498 !if(bband==kband) then 499 !err_ovlp=abs(mag_ovlp-1.0) 500 !else 501 !err_ovlp=abs(mag_ovlp) 502 !end if 503 !max_err_ovlp=MAX(max_err_ovlp,err_ovlp) 504 !end do ! end loop over kband 505 !end do ! end loop over bband 506 !write(std_out,'(a,i4,es16.8)')' JWZ Debug: berrphase_new ikpt ovlp err : ',& 507 !& ikpt,max_err_ovlp 508 !ABI_FREE(bra) 509 !ABI_FREE(ket) 510 ! 511 !!========================================= 512 !! end code to test orthonormality of cg_k 513 !!========================================= 514 515 pel(:) = zero ; pelev(:)=zero ; pion(:) = zero ; ptot(:)=zero ; red_ptot(:)=zero 516 517 minbd = 1 ; maxbd = dtefield%mband_occ 518 519 if(calc_epaw3_force) dtefield%epawf3(:,:,:) = zero 520 if(calc_epaw3_stress) dtefield%epaws3(:,:,:) = zero 521 522 do idir = 1, 3 523 524 ! dtm(:,:) = zero 525 dtm_mult(:,:,:) = zero 526 if (calc_epaw3_force) dsdr_sum(:,:,:) = zero 527 if (calc_epaw3_stress) dsds_sum(:,:,:) = zero 528 529 if (dtset%rfdir(idir) /= 1) cycle 530 531 if (abs(dtefield%efield_dot(idir)) < tol12) dtefield%sflag(:,:,:,idir) = 0 532 533 ! calculate vector steps in k space 534 dk(:) = dtefield%dkvecs(:,idir) 535 gpard(:) = dk(1)*gprimd(:,1) + dk(2)*gprimd(:,2) + dk(3)*gprimd(:,3) 536 gmod = sqrt(dot_product(gpard,gpard)) 537 538 write(message,'(a,a,a,3f9.5,a,a,3f9.5,a)')ch10,& 539 & ' Computing the polarization (Berry phase) for reciprocal vector:',ch10,& 540 & dk(:),' (in reduced coordinates)',ch10,& 541 & gpard(1:3),' (in cartesian coordinates - atomic units)' 542 call wrtout(std_out,message,'COLL') 543 if (unit_out /= 0) then 544 call wrtout(unit_out,message,'COLL') 545 end if 546 547 write(message,'(a,i5,a,a,i5)')& 548 & ' Number of strings: ',dtefield%nstr(idir),ch10,& 549 & ' Number of k points in string:', dtefield%nkstr(idir) 550 call wrtout(std_out,message,'COLL') 551 if (unit_out /= 0) then 552 call wrtout(unit_out,message,'COLL') 553 end if 554 555 ! Check whether the polarization or the ddk must be computed 556 557 ! nfor = 1 : to compute P, I only need the WF at k + dk 558 ! nfor = 2 : to compute the ddk, I need the WF at k + dk and k - dk 559 ! dkinv : +-1/2dk 560 561 562 ! default for polarization 563 nfor = 1 564 if (ddkflag == 1) then 565 nfor = 2 566 end if 567 568 if (ddkflag == 1) then 569 570 cg1(:,:) = zero 571 dkinv = one/(two*dk(idir)) 572 573 write(message,'(a,a,a,3f9.5,a,a,3f9.5,a)')ch10,& 574 & ' Computing the ddk (Berry phase) for reciprocal vector:',ch10,& 575 & dk(:),' (in reduced coordinates)',ch10,& 576 & gpard(1:3),' (in cartesian coordinates - atomic units)' 577 call wrtout(std_out,message,'COLL') 578 if (unit_out /= 0) then 579 call wrtout(unit_out,message,'COLL') 580 end if 581 582 end if 583 584 ! From smatrix routine: det_inv_smat = type of calculation 585 ! 1 : compute inverse of the overlap matrix 586 ! 10 : compute determinant of the overlap matrix 587 ! 11 : compute determinant and inverse of the overlap matrix 588 inv_smat = 0 589 det_smat = 0 590 591 ! for ddk need inverse matrix 592 if (ddkflag == 1) then 593 inv_smat = 1 594 end if 595 596 ! if polarization is requested need smat determinant as well 597 if (polflag == 1) then 598 det_smat = 1 599 end if 600 601 ! electric fields with PAW also needs S_inverse for forces and stresses, even just for polarization 602 if (calc_epaw3_force .or. calc_epaw3_stress) then 603 inv_smat = 1 604 end if 605 606 det_inv_smat = 10*det_smat + inv_smat 607 608 !-------------------------------------------------------------------- 609 ! for each dk we require, calculate the smatrix, derivatives etc... 610 !-------------------------------------------------------------------- 611 do ifor = 1, nfor 612 613 if (ifor == 2) then 614 dk(:) = -1_dp*dk(:) 615 ! only the inverse of the overlap matrix is required on second pass, speeds things up a bit 616 det_inv_smat = 1 617 dkinv = -1_dp*dkinv 618 end if 619 620 621 ! Compute the determinant and/or the inverse of the overlap matrix 622 ! for each pair of k-points < u_nk | u_nk+dk > 623 624 icg = 0 ; icg1 = 0 625 ABI_MALLOC(smat_k,(2,dtefield%mband_occ,dtefield%mband_occ)) 626 ABI_MALLOC(smat_inv,(2,dtefield%mband_occ,dtefield%mband_occ)) 627 ABI_MALLOC(smat_k_paw,(2,usepaw*dtefield%mband_occ,usepaw*dtefield%mband_occ)) 628 if (calc_epaw3_force .or. calc_epaw3_stress) then ! dsdr needed for forces and stresses in electric field with PAW 629 ABI_MALLOC(dsdr,(2,natom,ncpgr,usepaw*dtefield%mband_occ,usepaw*dtefield%mband_occ)) 630 dsdr = zero 631 end if 632 633 634 ! Loop on the values of ikpt_loc and ikpt1 : 635 ! ikpt1 is incremented one by one, and number the k points in the FBZ 636 ! ikpt1i refer to the k point numbering in the IBZ 637 ! ikpt_loc differs from ikpt1 only in the parallel case, and gives 638 ! the index of the k point in the FBZ, in the set treated by the present processor 639 ! NOTE : in order to allow synchronisation, ikpt_loc contain information about 640 ! ikpt AND ISPPOL ! 641 ! It means that the following loop is equivalent to a double loop : 642 ! do isppol = 1, nsppol 643 ! do ikpt1 = 1, dtefield%fmkmem 644 ! 645 do ikpt_loc = 1, dtefield%fmkmem_max*nsppol 646 647 ikpt1=mpi_enreg%kpt_loc2fbz_sp(me, ikpt_loc,1) 648 isppol=mpi_enreg%kpt_loc2fbz_sp(me, ikpt_loc,2) 649 650 ! if this k and spin are for me do it 651 if (ikpt1 > 0 .and. isppol > 0) then 652 653 ikpt1i = dtefield%indkk_f2ibz(ikpt1,1) 654 nband_k = dtset%nband(ikpt1i + (isppol-1)*dtset%nkpt) 655 656 ! DEBUG 657 ! Please keep this debugging feature 658 ! write(std_out,'(a,5i4)' )' berryphase_new : ikpt_loc,ikpt1,isppol,idir,ifor=',& 659 ! & ikpt_loc,ikpt1,isppol,idir,ifor 660 ! ENDDEBUG 661 662 inibz=0 663 if (dtset%kptns(1,ikpt1i) == dtefield%fkptns(1,ikpt1) .and. & 664 & dtset%kptns(2,ikpt1i) == dtefield%fkptns(2,ikpt1) .and. & 665 & dtset%kptns(3,ikpt1i) == dtefield%fkptns(3,ikpt1)) inibz=1 666 667 ikg = dtefield%fkgindex(ikpt1) 668 ! ikpt2 = dtefield%ikpt_dk(ikpt1,ifor,idir) 669 ! ikpt2i = dtefield%indkk_f2ibz(ikpt2,1) 670 671 ! ikpt3(istep) : index of kpt1 + istep*dk in the FBZ 672 ! ikpt3i(istep) : index of kpt1 + istep*dk in the IBZ 673 ikpt3(1) = dtefield%ikpt_dk(ikpt1,ifor,idir) 674 ikpt3i(1) = dtefield%indkk_f2ibz(ikpt3(1),1) 675 do istep = 1, berrystep-1 676 ikpt3(istep+1) = dtefield%ikpt_dk(ikpt3(istep),ifor,idir) 677 ikpt3i(istep+1) = dtefield%indkk_f2ibz(ikpt3(istep+1),1) 678 end do 679 680 ! itrs = 0 681 ! if (dtefield%indkk_f2ibz(ikpt1,6) == 1 ) itrs = itrs + 1 682 ! if (dtefield%indkk_f2ibz(ikpt2,6) == 1 ) itrs = itrs + 10 683 684 itrs_mult(:)=0 685 if (dtefield%indkk_f2ibz(ikpt1,6) == 1 ) itrs_mult(:) = itrs_mult(:) + 1 686 do istep=1,berrystep 687 if (dtefield%indkk_f2ibz(ikpt3(istep),6) == 1 ) itrs_mult(istep) = itrs_mult(istep) + 10 688 end do 689 690 npw_k1 = npwarr(ikpt1i) 691 ! npw_k2 = npwarr(ikpt2i) 692 693 do istep = 1, berrystep 694 npw_k3(istep)=npwarr(ikpt3i(istep)) 695 end do 696 697 ! ji: the loop is over the FBZ, but sflag and smat only apply to the IBZ 698 if ( efield_flag .and. inibz == 1) then !!HONG 699 ikpt1i_sp=ikpt1i+(isppol-1)*dtset%nkpt 700 smat_k(:,:,:) = dtefield%smat(:,:,:,ikpt1i_sp,ifor,idir) 701 else 702 smat_k(:,:,:) = zero 703 end if 704 705 ! pwind_k(1:npw_k1) = pwind(ikg+1:ikg+npw_k1,ifor,idir) 706 pwnsfac_k(1,1:npw_k1) = pwnsfac(1,ikg+1:ikg+npw_k1) 707 pwnsfac_k(2,1:npw_k1) = pwnsfac(2,ikg+1:ikg+npw_k1) 708 709 ! the array needed to compute the overlap matrix between k and k+istep*dk (with multiple steps) 710 ! the 0-case (no corresponding pw in k and k+dk) could be handled better (k+2*dk could have a corresponding pw ?) 711 pwind_k_mult(1:npw_k1,1)=pwind(ikg+1:ikg+npw_k1,ifor,idir) 712 do istep = 1, berrystep-1 713 do jj=1, npw_k1 714 if(pwind_k_mult(jj,istep)/=0)then 715 pwind_k_mult(jj,istep+1) = pwind(dtefield%fkgindex(ikpt3(istep))+pwind_k_mult(jj,istep),ifor,idir) 716 else 717 pwind_k_mult(jj,istep+1) = 0 718 end if 719 end do 720 end do 721 722 ! DEBUG 723 ! write(std_out,*)' berryphase_new : dtset%berryopt,inibz,ikpt1i,isppol,dtset%nkpt,ifor,idir', & 724 ! & dtset%berryopt,inibz,ikpt1i,isppol,dtset%nkpt,ifor,idir 725 ! write(std_out,'(a,4i4)' )' berryphase_new : sflag_k(:)=',sflag_k(:) 726 ! ENDDEBUG 727 728 if ( efield_flag .and. inibz == 1) then !!HONG 729 ikpt1i_sp=ikpt1i+(isppol-1)*dtset%nkpt 730 sflag_k(:) = dtefield%sflag(:,ikpt1i_sp,ifor,idir) 731 else 732 sflag_k(:) = 0 733 end if 734 735 if (usepaw == 1) then 736 icp1=dtefield%cprjindex(ikpt1i,isppol) 737 call pawcprj_get(atindx1,cprj_k,cprj,natom,1,icp1,ikpt1i,0,isppol,& 738 & mband,mkmem,natom,dtefield%mband_occ,dtefield%mband_occ,& 739 & my_nspinor,nsppol,0,mpicomm=mpi_enreg%comm_kpt,& 740 & proc_distrb=mpi_enreg%proc_distrb) 741 742 if ( ikpt1i /= ikpt1 ) then 743 call pawcprj_copy(cprj_k,cprj_ikn) 744 call pawcprj_symkn(cprj_fkn,cprj_ikn,dtefield%atom_indsym,dimlmn,-1,indlmn,& 745 & dtefield%indkk_f2ibz(ikpt1,2),dtefield%indkk_f2ibz(ikpt1,6),& 746 & dtefield%fkptns(:,dtefield%i2fbz(ikpt1i)),& 747 & dtefield%lmax,dtefield%lmnmax,mband,natom,dtefield%mband_occ,my_nspinor,& 748 & dtefield%nsym,ntypat,typat,dtefield%zarot) 749 call pawcprj_copy(cprj_fkn,cprj_k) 750 end if 751 752 end if ! end if usepaw 753 754 ! DEBUG 755 ! write(std_out,'(a,4i4)' )' berryphase_new : sflag_k(:)=',sflag_k(:) 756 ! ENDDEBUG 757 758 ! DEBUG 759 ! write(std_out,'(a,7i4)')'me, idir,ifor, ikpt_loc, ikpt1, isppol = ',& 760 ! & me,idir,ifor,ikpt_loc,ikpt1,isppol 761 ! write(std_out,'(a,10i3)')'pwind_k(1:10) = ',pwind_k(1:10) 762 ! ENDDEBUG 763 764 do istep=1,berrystep 765 sflag_k_mult(:,istep) = sflag_k(:) 766 end do 767 768 end if ! end check that ikpt1 > 0 and isppol > 0 769 770 ! -------------------------------------------------------------------------------- 771 ! Communication 772 ! -------------------------------------------------------------------------------- 773 774 do istep=1,berrystep 775 776 ! if(ikpt_loc <= nsppol*dtefield%fmkmem) then 777 if (ikpt1 > 0 .and. isppol > 0) then ! I currently have a true kpt to use 778 779 count = npw_k3(istep)*my_nspinor*nband_k 780 ABI_MALLOC(cgq,(2,count)) 781 cgq = zero 782 source = me 783 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt3i(istep),1,nband_k,isppol,me)) then 784 ! I need the datas from someone else 785 source = mpi_enreg%proc_distrb(ikpt3i(istep),1,isppol) 786 end if 787 else 788 source = -1 ! I do not have a kpt to use 789 end if 790 791 do dest = 0, nproc-1 792 793 if ((dest==me) .and. (ikpt1>0) .and. (isppol>0)) then 794 ! I am destination and I have something to do 795 ! if (mpi_enreg%paral_compil_kpt == 1) write(std_out,*) & 796 ! & 'coucou 2, mpi_enreg%proc_distrb(ikpt3i(istep),1:nband_k,isppol) : ', & 797 ! & mpi_enreg%proc_distrb(ikpt3i(istep),1:nband_k,isppol) 798 ! write(std_out,*)'ikpt3i(istep) ', ikpt3i(istep) 799 ! write(std_out,*)'nband_k ',nband_k 800 ! write(std_out,*)'isppol ', isppol 801 ! write(std_out,*)'mpi_enreg%proc_distrb',mpi_enreg%proc_distrb 802 803 if (source == me) then 804 ! I am destination and source 805 ! DEBUG 806 ! write(std_out,*)'copying ... ' 807 ! write(std_out,*)'me: ',me, 'ikpt3i(istep) ', ikpt3i(istep), 'isppol ', isppol 808 ! ENDDEBUG 809 810 ! pwnsfac 811 idum = dtefield%fkgindex(ikpt3(istep)) 812 pwnsfac_k(3,1:npw_k3(istep)) = pwnsfac(1,idum+1:idum+npw_k3(istep)) 813 pwnsfac_k(4,1:npw_k3(istep)) = pwnsfac(2,idum+1:idum+npw_k3(istep)) 814 815 ! cgq (and cprj) 816 icg1 = dtefield%cgindex(ikpt3i(istep),isppol) 817 818 if (usepaw == 1) then 819 icp2=dtefield%cprjindex(ikpt3i(istep),isppol) 820 call pawcprj_get(atindx1,cprj_kb,cprj,natom,1,icp2,ikpt3i(istep),0,isppol,& 821 & mband,mkmem,natom,dtefield%mband_occ,dtefield%mband_occ,my_nspinor,& 822 & nsppol,0,mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 823 end if 824 825 cgq(:,1:count) = cg(:,icg1+1:icg1+count) 826 ! if (usepaw == 1) then 827 ! call pawcprj_copy(cprj_buf,cprj_kb) 828 ! end if 829 830 ! if ((source /= me)) then 831 else 832 ! I am the destination but not the source -> receive 833 ! DEBUG 834 ! write(std_out,'(a)')'receiving ...' 835 ! write(std_out,'(a,i4,a,i4,a,i4,a,i4)')'me: ',me, 'source ', source,'ikpt3i(istep) ', ikpt3i(istep), 'isppol ', isppol 836 ! ENDDEBUG 837 838 ! receive pwnsfac 839 ABI_MALLOC(buffer,(2,npw_k3(istep))) 840 tag = ikpt3(istep) + (isppol - 1)*dtefield%fnkpt 841 call xmpi_recv(buffer,source,tag,spaceComm,ierr) 842 pwnsfac_k(3,1:npw_k3(istep)) = buffer(1,1:npw_k3(istep)) 843 pwnsfac_k(4,1:npw_k3(istep)) = buffer(2,1:npw_k3(istep)) 844 ABI_FREE(buffer) 845 846 ! receive cgq (and cprj) 847 tag = ikpt3i(istep) + (isppol - 1)*nkpt 848 call xmpi_recv(cgq,source,tag,spaceComm,ierr) 849 850 if (usepaw == 1) then 851 call pawcprj_mpi_recv(natom,n2dim,dimlmn,ncpgr,cprj_kb,source,spaceComm,ierr) 852 end if 853 854 end if 855 856 else if (dest /= me) then 857 858 ! jkpt is the kpt which is being treated by dest 859 ! jsppol is his isppol 860 jkpt = mpi_enreg%kpt_loc2fbz_sp(dest, ikpt_loc,1) 861 jsppol = mpi_enreg%kpt_loc2fbz_sp(dest, ikpt_loc,2) 862 863 if (jkpt > 0 .and. jsppol > 0) then ! dest is treating a true kpt 864 865 jkpt2 = dtefield%ikpt_dk(jkpt,ifor,idir) 866 jkpt2i = dtefield%indkk_f2ibz(jkpt2,1) 867 868 ! check if I am his source 869 if((mpi_enreg%proc_distrb(jkpt2i,1,jsppol) == me)) then 870 ! I know something about jkpt3i and I must send it 871 ! DEBUG 872 ! write(std_out,'(a)')'sending ...' 873 ! write(std_out,'(a,i4,a,i4,a,i4,a,i4)')'dest: ',dest,' me: ',me,& 874 ! & ' jkpt2i ',jkpt2i,' jsppol: ',jsppol 875 ! ENDDEBUG 876 877 ! pwnsfac 878 tag = jkpt2 + (jsppol - 1)*dtefield%fnkpt 879 count1 = npwarr(jkpt2i) 880 ABI_MALLOC(buffer,(2,count1)) 881 idum = dtefield%fkgindex(jkpt2) 882 buffer(1,1:count1) = pwnsfac(1,idum+1:idum+count1) 883 buffer(2,1:count1) = pwnsfac(2,idum+1:idum+count1) 884 call xmpi_send(buffer,dest,tag,spaceComm,ierr) 885 ABI_FREE(buffer) 886 887 ! cgq (and cprj) 888 icg1 = dtefield%cgindex(jkpt2i,jsppol) 889 890 if (usepaw == 1) then 891 icp2=dtefield%cprjindex(jkpt2i,jsppol) 892 call pawcprj_get(atindx1,cprj_buf,cprj,natom,1,icp2,jkpt2i,0,jsppol,& 893 & mband,mkmem,natom,dtefield%mband_occ,dtefield%mband_occ,& 894 & my_nspinor,nsppol,0,mpicomm=mpi_enreg%comm_kpt,& 895 & proc_distrb=mpi_enreg%proc_distrb) 896 end if 897 898 tag = jkpt2i + (jsppol - 1)*nkpt 899 count1 = npwarr(jkpt2i)*my_nspinor*nband_k 900 ABI_MALLOC(buffer,(2,count1)) 901 buffer(:,1:count1) = cg(:,icg1+1:icg1+count1) 902 call xmpi_send(buffer,dest,tag,spaceComm,ierr) 903 ABI_FREE(buffer) 904 905 if (usepaw == 1 ) then 906 call pawcprj_mpi_send(natom,n2dim,dimlmn,ncpgr,cprj_buf,dest,spaceComm,ierr) 907 end if 908 909 end if ! end check that I am his source 910 911 end if ! end check that jkpt > 0 and jsppol > 0 912 913 end if ! end if statements on dest == me or dest /= me 914 915 end do ! end loop over dest = 0, nproc - 1 916 917 if (ikpt1 > 0 .and. isppol > 0) then ! if I am treating a kpt, compute the smatrix 918 919 if (usepaw == 1) then 920 if (ikpt3(istep) /= ikpt3i(istep)) then ! cprj_kb refers to ikpt3i(istep), must compute ikpt3(istep) value 921 call pawcprj_copy(cprj_kb,cprj_ikn) 922 923 call pawcprj_symkn(cprj_fkn,cprj_ikn,dtefield%atom_indsym,dimlmn,-1,indlmn,& 924 & dtefield%indkk_f2ibz(ikpt3(istep),2),dtefield%indkk_f2ibz(ikpt3(istep),6),& 925 & dtefield%fkptns(:,dtefield%i2fbz(ikpt3i(istep))),& 926 & dtefield%lmax,dtefield%lmnmax,mband,natom,& 927 & dtefield%mband_occ,my_nspinor,dtefield%nsym,ntypat,typat,& 928 & dtefield%zarot) 929 call pawcprj_copy(cprj_fkn,cprj_kb) 930 end if 931 call smatrix_k_paw(cprj_k,cprj_kb,dtefield,idir,ifor,mband,natom,smat_k_paw,typat) 932 ! write(std_out,'(a,5i4)')' JWZ berryphase_new : ikpt_loc,ikpt1,ikpt1i,ikpt2,ikpt2i ',ikpt_loc,ikpt1,ikpt1i,ikpt3(istep),ikpt3i(istep) 933 ! call smatrix_k0_paw(atindx1,cprj_k,cprj_k,dtefield,ikpt1i,idir,ifor,& 934 ! & mband,mpi_enreg,natom,ntypat,pawtab,smat_k_paw,typat) 935 if (calc_epaw3_force .or. calc_epaw3_stress) then 936 call dsdr_k_paw(cprj_k,cprj_kb,dsdr,dtefield,idir,ifor,mband,natom,ncpgr,typat) 937 end if 938 end if 939 940 icg1 = 0 941 icg = dtefield%cgindex(ikpt1i,isppol) 942 ! DEBUG 943 ! if(istep<=2)then 944 ! if(ikpt1==1)then 945 ! write(std_out,'(a,2i4,3e15.4)')'istep ikpt3, kpt, cgq', istep, ikpt3(istep), dtefield%fkptns(:,ikpt3(istep)) 946 ! write(std_out,*) cgq 947 ! write(std_out,*) 948 ! end if 949 ! end if 950 ! ENDDEBUG 951 call smatrix(cg,cgq,cg1_k,ddkflag,dtm_k,icg,icg1,itrs_mult(istep),det_inv_smat,maxbd,& 952 & mcg,count,mcg1_k,minbd,& 953 & mpw,dtefield%mband_occ,dtefield%nband_occ(isppol),& 954 & npw_k1,npw_k3(istep),my_nspinor,pwind_k_mult(:,istep),pwnsfac_k,sflag_k_mult(:,istep),& 955 & shiftbd,smat_inv,smat_k,smat_k_paw,usepaw) 956 957 ! in finite electric field case with paw must save additional F3 term in forces 958 if(calc_epaw3_force) then 959 ! when ncpgr = 3, gradients are wrt to atom displacements 960 ! but when ncpgr = 9, first 6 gradients are wrt strains, last three are displacements 961 icpgr_offset = 0 962 if (ncpgr == 9) icpgr_offset = 6 963 do iatom = 1, natom 964 do fdir = 1, 3 965 dsdr_sum(iatom,fdir,ikpt1+(isppol-1)*dtefield%fnkpt) = zero 966 do iband = 1, dtefield%nband_occ(isppol) 967 do jband = 1, dtefield%nband_occ(isppol) 968 ! collect Im{Trace{S^{-1}.dS/dR}} for this k point 969 dsdr_sum(iatom,fdir,ikpt1+(isppol-1)*dtefield%fnkpt) = & 970 & dsdr_sum(iatom,fdir,ikpt1+(isppol-1)*dtefield%fnkpt) + & 971 & smat_inv(2,iband,jband)*dsdr(1,iatom,icpgr_offset+fdir,jband,iband) + & 972 & smat_inv(1,iband,jband)*dsdr(2,iatom,icpgr_offset+fdir,jband,iband) 973 end do ! end sum over jband 974 end do ! end sum over iband 975 end do ! end sum over fdir 976 end do ! end sum over iatom 977 end if ! end check on calc_epaw3_force 978 979 ! in finite electric field case with paw must save additional F3 term in stress 980 ! note that when strains are present they are always saved before forces 981 ! therefore no need for icpgr_offset in this case 982 if(calc_epaw3_stress) then 983 do iatom = 1, natom 984 do fdir = 1, 6 985 dsds_sum(iatom,fdir,ikpt1+(isppol-1)*dtefield%fnkpt) = zero 986 do iband = 1, dtefield%nband_occ(isppol) 987 do jband = 1, dtefield%nband_occ(isppol) 988 ! collect Im{Trace{S^{-1}.dS/de}} for this k point 989 dsds_sum(iatom,fdir,ikpt1+(isppol-1)*dtefield%fnkpt) = & 990 & dsds_sum(iatom,fdir,ikpt1+(isppol-1)*dtefield%fnkpt) + & 991 & smat_inv(2,iband,jband)*dsdr(1,iatom,fdir,jband,iband) + & 992 & smat_inv(1,iband,jband)*dsdr(2,iatom,fdir,jband,iband) 993 end do ! end sum over jband 994 end do ! end sum over iband 995 end do ! end sum over fdir 996 end do ! end sum over iatom 997 end if ! end check on calc_epaw3_stress 998 999 if ((det_inv_smat == 10).or.(det_inv_smat == 11)) then 1000 1001 if (sqrt(dtm_k(1)*dtm_k(1) + dtm_k(2)*dtm_k(2)) < tol12) then 1002 ! EB: the ABI_BUG has been replaced here by what is done in 67_common/m_cgwf.F90 1003 ! This avoid the code to stop for phonons under E-field too. 1004 ! TODO: Since the same is done in m_cgwf.F90 and in m_berryphase_new.F90, 1005 ! rationalization should be done with one single module. 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. Fixing...' 1009 ! Try this: 1010 write(std_out,*)message,dtm_k(1:2) 1011 if(abs(dtm_k(1))<=1d-12)dtm_k(1)=1d-12 1012 if(abs(dtm_k(2))<=1d-12)dtm_k(2)=1d-12 1013 write(std_out,*)' Changing to:',dtm_k(1:2) 1014 ! ABI_BUG(message) 1015 end if 1016 1017 dtm_mult(1,ikpt1+(isppol-1)*dtefield%fnkpt,istep) = dtm_k(1) 1018 dtm_mult(2,ikpt1+(isppol-1)*dtefield%fnkpt,istep) = dtm_k(2) 1019 1020 end if 1021 1022 if ( efield_flag .and. inibz == 1 .and. istep == 1) then !!HONG 1023 ikpt1i_sp=ikpt1i+(isppol-1)*dtset%nkpt 1024 dtefield%smat(:,:,:,ikpt1i_sp,ifor,idir) = & 1025 & smat_k(:,:,:) 1026 dtefield%sflag(:,ikpt1i_sp,ifor,idir) = & 1027 & sflag_k_mult(:,1) 1028 end if 1029 1030 ! for IBZ k-points and first step, add 1031 if ((ddkflag==1 .and.((det_inv_smat == 1).or.(det_inv_smat == 11))) .and. inibz == 1 .and. istep == 1) then 1032 cg1(:,icg + 1: icg + npw_k1*my_nspinor*nband_k) = & 1033 cg1(:,icg + 1:icg + npw_k1*my_nspinor*nband_k) + & 1034 dkinv*cg1_k(:,1:npw_k1*my_nspinor*nband_k) 1035 end if 1036 1037 ABI_FREE(cgq) 1038 1039 end if ! end if ikpt1 > 0 and isppol > 0 1040 1041 end do ! end loop over istep 1042 1043 ! if (ikpt_loc <= dtefield%fmkmem) sflag_k(:) = sflag_k_mult(:,1) 1044 if (ikpt1 > 0) sflag_k(:) = sflag_k_mult(:,1) 1045 1046 end do ! close loop over ikpt_loc (k-points, isppol) 1047 1048 ABI_FREE(smat_inv) 1049 ABI_FREE(smat_k) 1050 ABI_FREE(smat_k_paw) 1051 if (calc_epaw3_force .or. calc_epaw3_stress) then 1052 ABI_FREE(dsdr) 1053 end if 1054 1055 end do ! close loop over ifor 1056 1057 ! MPI communicate stuff between everyone 1058 if (nproc>1) then 1059 count = 2*dtefield%fnkpt*nsppol*berrystep 1060 ABI_MALLOC(buffer1,(count)) 1061 ABI_MALLOC(buffer2,(count)) 1062 buffer1(:) = reshape(dtm_mult(:,:,:),(/count/)) 1063 call xmpi_sum(buffer1,buffer2,count,spaceComm,ierr) 1064 dtm_mult(:,:,:) = reshape(buffer2(:),(/2,dtefield%fnkpt*nsppol,berrystep/)) 1065 ABI_FREE(buffer1) 1066 ABI_FREE(buffer2) 1067 if (calc_epaw3_force) then 1068 count = natom*3*dtefield%fnkpt*nsppol 1069 ABI_MALLOC(buffer1,(count)) 1070 ABI_MALLOC(buffer2,(count)) 1071 buffer1(:) = reshape(dsdr_sum(:,:,:),(/count/)) 1072 call xmpi_sum(buffer1,buffer2,count,spaceComm,ierr) 1073 dsdr_sum(:,:,:) = reshape(buffer2(:),(/natom,3,dtefield%fnkpt*nsppol/)) 1074 ABI_FREE(buffer1) 1075 ABI_FREE(buffer2) 1076 end if 1077 if (calc_epaw3_stress) then 1078 count = natom*6*dtefield%fnkpt*nsppol 1079 ABI_MALLOC(buffer1,(count)) 1080 ABI_MALLOC(buffer2,(count)) 1081 buffer1(:) = reshape(dsds_sum(:,:,:),(/count/)) 1082 call xmpi_sum(buffer1,buffer2,count,spaceComm,ierr) 1083 dsds_sum(:,:,:) = reshape(buffer2(:),(/natom,6,dtefield%fnkpt*nsppol/)) 1084 ABI_FREE(buffer1) 1085 ABI_FREE(buffer2) 1086 end if 1087 end if ! if parallel 1088 1089 ! DEBUG 1090 ! write(std_out,*) 1091 ! write(std_out,*)'istep = 1, nsppol =',nsppol 1092 ! istep=1 1093 ! isppol=1 1094 ! do jkpt = 1, dtefield%fnkpt 1095 ! 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) 1096 ! end do 1097 ! write(std_out,*) 1098 ! write(std_out,*) "istep = 2" 1099 ! if(berrystep>=2)then 1100 ! istep=2 1101 ! isppol=1 1102 ! do jkpt = 1, dtefield%fnkpt 1103 ! 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) 1104 ! end do 1105 ! end if 1106 ! ENDDEBUG 1107 1108 ! =========================================================================== 1109 ! in DDK case everything has been calculated above from finite difference 1110 ! Now write the ddk WF to a file 1111 ! =========================================================================== 1112 1113 if (ddkflag == 1) then 1114 1115 pertcase = idir + 3*natom 1116 response = 1 1117 call appdig(pertcase,dtfil%fnameabo_1wf,fiwf1o) 1118 ABI_MALLOC(resid,(mband*nkpt*nsppol)) 1119 resid(:) = zero 1120 1121 call outresid(dtset,dtset%kptns,mband,& 1122 & dtset%nband,nkpt,& 1123 & nsppol,resid) 1124 1125 call outwf(cg1,dtset,psps,eig_dum,fiwf1o,hdr,kg,dtset%kptns,& 1126 & mband,mcg,mkmem,mpi_enreg,mpw,natom,dtset%nband,& 1127 & nkpt,npwarr,nsppol,& 1128 & occ_dum,response,dtfil%unwff2,wfs,wvl) 1129 1130 if (save_cg13) then 1131 cg13(1:2,1:mcg13,idir) = cg1(1:2,1:mcg) 1132 end if 1133 1134 ABI_FREE(resid) 1135 end if ! ddkflag == 1 1136 ! end of ddk part for this idir 1137 1138 1139 ! =========================================================================== 1140 ! Compute the Berry phase polarization 1141 ! =========================================================================== 1142 1143 if (polflag == 1) then 1144 1145 ! Compute the electronic Berry phase 1146 1147 polb_mult(:,:)=zero 1148 do istep = 1,berrystep 1149 1150 if(berrystep==1) then 1151 write(message,'(a,a)')ch10,& 1152 & ' Compute the electronic contribution to polarization' 1153 call wrtout(std_out,message,'COLL') 1154 else 1155 write(message,'(a,a,i4,a)')ch10,& 1156 & ' Compute the electronic contribution to polarization for a step of istep=',& 1157 & istep,'*dk' 1158 call wrtout(std_out,message,'COLL') 1159 end if 1160 1161 if(istep /= 1) then 1162 ! construct the strings for a step of istep*dk 1163 ! string length 1164 istr=1 1165 nkstr=1 1166 ikpt1=1 1167 do ikpt=1,dtefield%fnkpt 1168 do jstep = 1,istep 1169 ikpt1 = dtefield%ikpt_dk(ikpt1,1,idir) 1170 end do 1171 if (ikpt1 == 1) exit 1172 nkstr = nkstr + 1 1173 end do 1174 ! Check that the string length is a divisor of nkpt 1175 if(mod(dtefield%fnkpt,nkstr) /= 0) then 1176 write(message,'(a,i5,a,i5,a,i7)')& 1177 & ' For istep = ', istep,& 1178 & ' The string length = ',nkstr,& 1179 & ', is not a divisor of fnkpt =',dtefield%fnkpt 1180 ABI_BUG(message) 1181 end if 1182 nstr = dtefield%fnkpt/nkstr 1183 1184 write(message,'(a,i1,a,i2,a,i3,a,i6)')& 1185 & ' berryphase_new: for direction ',idir, ' and istep ', istep, ', nkstr = ',nkstr,& 1186 & ', nstr = ',nstr 1187 call wrtout(std_out,message,'COLL') 1188 call wrtout(ab_out,message,'COLL') 1189 1190 ABI_MALLOC(idxkstr_mult,(nkstr,nstr)) 1191 iunmark = 1 1192 kpt_mark(:)=0 1193 do istr=1,nstr 1194 do while(kpt_mark(iunmark) /= 0) 1195 iunmark = iunmark + 1 1196 end do 1197 idxkstr_mult(1,istr) = iunmark 1198 kpt_mark(iunmark)=1 1199 1200 ikpt1 = idxkstr_mult(1,istr) 1201 do jkstr=2, nkstr 1202 do jstep = 1, istep 1203 ikpt1 = dtefield%ikpt_dk(ikpt1,1,idir) 1204 end do 1205 idxkstr_mult(jkstr,istr) = ikpt1 1206 kpt_mark(ikpt1) = 1 1207 end do 1208 end do 1209 else 1210 nstr = dtefield%nstr(idir) 1211 nkstr = dtefield%nkstr(idir) 1212 ABI_MALLOC(idxkstr_mult,(nkstr,nstr)) 1213 idxkstr_mult(:,:) = dtefield%idxkstr(1:nkstr,1:nstr,idir) 1214 end if 1215 ! DEBUG 1216 ! do istr=1,nstr 1217 ! write(std_out,*)'string ', idxkstr_mult(:,istr) 1218 ! end do 1219 ! ENDBEBUG 1220 1221 ABI_MALLOC(det_string,(2,nstr)) 1222 ABI_MALLOC(polberry,(nstr)) 1223 write(message,'(a,10x,a,10x,a)')ch10,& 1224 & 'istr','polberry(istr)' 1225 call wrtout(std_out,message,'COLL') 1226 1227 polbtot = zero 1228 do isppol = 1, nsppol 1229 1230 det_string(1,:) = one ; det_string(2,:) = zero 1231 dtm_k(:) = one 1232 det_average(:) = zero 1233 1234 1235 do istr = 1, nstr 1236 1237 if(calc_epaw3_force) epawf3_str(:,:,:) = zero 1238 if(calc_epaw3_stress) epaws3_str(:,:,:) = zero 1239 1240 do jkstr = 1, nkstr 1241 1242 ikpt=idxkstr_mult(jkstr,istr) 1243 1244 dtm_real=dtm_mult(1,ikpt+(isppol-1)*dtefield%fnkpt,istep) 1245 dtm_imag=dtm_mult(2,ikpt+(isppol-1)*dtefield%fnkpt,istep) 1246 1247 dtm_k(1) = det_string(1,istr)*dtm_real - & 1248 & det_string(2,istr)*dtm_imag 1249 dtm_k(2) = det_string(1,istr)*dtm_imag + & 1250 & det_string(2,istr)*dtm_real 1251 det_string(1:2,istr) = dtm_k(1:2) 1252 ! DEBUG 1253 ! write(std_out,'(a,i4,3e15.4,2e15.4)')'ikpt, kpt, dtm', ikpt, dtefield%fkptns(:,ikpt), dtm_k 1254 ! ENDDEBUG 1255 1256 if(calc_epaw3_force) then 1257 do iatom = 1, natom 1258 do fdir = 1, 3 1259 epawf3_str(iatom,idir,fdir) = epawf3_str(iatom,idir,fdir) + & 1260 & dsdr_sum(iatom,fdir,ikpt+(isppol-1)*dtefield%fnkpt) 1261 end do ! end loop over fdir 1262 end do ! end loop over natom 1263 end if ! end check on calc_epaw3_force 1264 if(calc_epaw3_stress) then 1265 do iatom = 1, natom 1266 do fdir = 1, 6 1267 epaws3_str(iatom,idir,fdir) = epaws3_str(iatom,idir,fdir) + & 1268 & dsds_sum(iatom,fdir,ikpt+(isppol-1)*dtefield%fnkpt) 1269 end do ! end loop over fdir 1270 end do ! end loop over natom 1271 end if ! end check on calc_epaw3_stress 1272 1273 end do 1274 1275 if(calc_epaw3_force) then 1276 do iatom = 1, natom 1277 do fdir = 1, 3 1278 dtefield%epawf3(iatom,idir,fdir) = dtefield%epawf3(iatom,idir,fdir) + & 1279 & epawf3_str(iatom,idir,fdir) 1280 end do ! end loop over fdir 1281 end do ! end loop over natom 1282 end if ! end check on calc_epaw3_force 1283 if(calc_epaw3_stress) then 1284 do iatom = 1, natom 1285 do fdir = 1, 6 1286 dtefield%epaws3(iatom,idir,fdir) = dtefield%epaws3(iatom,idir,fdir) + & 1287 & epaws3_str(iatom,idir,fdir) 1288 end do ! end loop over fdir 1289 end do ! end loop over natom 1290 end if ! end check on calc_epaw3_stress 1291 1292 det_average(:) = det_average(:) + & 1293 & det_string(:,istr)/dble(nstr) 1294 1295 end do 1296 1297 1298 ! correction to obtain a smooth logarithm of the determinant 1299 ABI_MALLOC(str_flag,(nstr)) 1300 ! DEBUG 1301 ! since we don't have any case of non-nul Chern number, 1302 ! we must change the det_string value "by brute force" if we want debug this 1303 ! allocate(det_string_test(2,dtefield%nstr(idir))) 1304 ! det_string_test(:,:)=det_string(:,:) 1305 ! kk=0 1306 ! det_string(1,1)=cos(2._dp*Pi*real(kk,dp)/four) 1307 ! det_string(2,1)=sin(2._dp*Pi*real(kk,dp)/four) 1308 ! jj=dtefield%str_neigh(1,1,idir) 1309 ! ll=dtefield%str_neigh(2,1,idir) 1310 ! do while (jj/=1) 1311 ! kk=kk+1 1312 ! det_string(1,jj)=cos(2._dp*Pi*real(kk,dp)/four) 1313 ! det_string(2,jj)=sin(2._dp*Pi*real(kk,dp)/four) 1314 ! det_string(1,ll)=cos(-2._dp*Pi*real(kk,dp)/four) 1315 ! det_string(2,ll)=sin(-2._dp*Pi*real(kk,dp)/four) 1316 ! jj=dtefield%str_neigh(1,jj,idir) 1317 ! ll=dtefield%str_neigh(2,ll,idir) 1318 ! enddo 1319 ! ENDDEBUG 1320 if (istep==1) then 1321 do ineigh_str = 1,2 1322 str_flag(:)=0 1323 delta_str(:) = & 1324 & dtefield%coord_str(:,dtefield%str_neigh(ineigh_str,1,idir),idir) - dtefield%coord_str(:,1,idir) 1325 dstr(:)= delta_str(:) - nint(delta_str(:)) - real(dtefield%strg_neigh(ineigh_str,1,:,idir),dp) 1326 dist_=0._dp 1327 do kk = 1,2 1328 do jj = 1,2 1329 dist_ = dist_ + dstr(kk)*dtefield%gmet_str(kk,jj,idir)*dstr(jj) 1330 end do 1331 end do 1332 dist_=sqrt(dist_) 1333 do istr = 1,dtefield%nstr(idir) 1334 if(str_flag(istr)==0)then 1335 ! write(std_out,*)'new string' 1336 str_flag(istr)=1 1337 call rhophi(det_string(:,istr),dphase,rho) 1338 ! write(std_out,'(i4,e15.4,e15.4,e15.4)')istr, det_string(:,istr),dphase 1339 dphase_init=dphase 1340 jstr = dtefield%str_neigh(ineigh_str,istr,idir) 1341 do while (istr/=jstr) 1342 str_flag(jstr)=1 1343 call rhophi(det_string(:,jstr),dphase_new,rho) 1344 jj=nint((dphase_new-dphase)/(2._dp*Pi)) 1345 ! DEBUG 1346 ! write(std_out,'(i4,e15.4,e15.4,e15.4,e15.4,i4)')jstr, det_string(:,jstr),dphase_new,dphase_new-dphase,jj 1347 ! ENDDEBUG 1348 dphase_new=dphase_new-two*Pi*real(jj,dp) 1349 if(jj/=0)then 1350 write(message,'(6a)') ch10,& 1351 & ' berryphase_new : WARNING -',ch10,& 1352 & ' the berry phase has some huge variation in the space of strings of k-points',ch10,& 1353 & ' ABINIT is trying to correct the berry phase, but it is highly experimental' 1354 call wrtout(std_out,message,'PERS') 1355 end if 1356 ! if(jj/=0)write(std_out,'(i4,e15.4,e15.4,e15.4,e15.4)')jstr, det_string(:,jstr),dphase_new,dphase_new-dphase 1357 dphase=dphase_new 1358 jstr=dtefield%str_neigh(ineigh_str,jstr,idir) 1359 end do 1360 ! write(std_out,*)dphase_init, dphase, (dphase-dphase_init)/(2._dp*Pi),nint((dphase-dphase_init)/(2._dp*Pi)) 1361 end if 1362 end do 1363 end do 1364 end if 1365 ABI_FREE(str_flag) 1366 ! DEBUG 1367 ! deallocate(dist_str) 1368 ! det_string(:,:)=det_string_test(:,:) 1369 ! deallocate(det_string_test) 1370 ! ENDDEBUG 1371 1372 ! First berry phase that corresponds to det_average 1373 ! phase0 = atan2(det_average(2),det_average(1)) 1374 call rhophi(det_average,phase0,rho) 1375 det_mod = det_average(1)**2+det_average(2)**2 1376 1377 ! Then berry phase that corresponds to each string relative to the average 1378 do istr = 1, nstr 1379 1380 rel_string(1) = (det_string(1,istr)*det_average(1) + & 1381 det_string(2,istr)*det_average(2))/det_mod 1382 rel_string(2) = (det_string(2,istr)*det_average(1) - & 1383 det_string(1,istr)*det_average(2))/det_mod 1384 ! dphase = atan2(rel_string(2),rel_string(1)) 1385 call rhophi(rel_string,dphase,rho) 1386 polberry(istr) = dtefield%sdeg*(phase0 + dphase)/two_pi 1387 polb_mult(isppol,istep) = polb_mult(isppol,istep) + polberry(istr)/(istep*dtefield%nstr(idir)) 1388 polb(isppol) = zero 1389 do jstep=1, istep 1390 polb(isppol)=polb(isppol)+coef(jstep,istep)*polb_mult(isppol,jstep) 1391 end do 1392 1393 write(message,'(10x,i6,7x,e16.9)')istr,polberry(istr) 1394 call wrtout(std_out,message,'COLL') 1395 1396 end do 1397 1398 if(berrystep>1)then 1399 write(message,'(9x,a,7x,e16.9,1x,a,i4,a,i4,a)')& 1400 & 'total',polb_mult(isppol,istep),'(isppol=',isppol,', istep=',istep,')'!,ch10 1401 call wrtout(std_out,message,'COLL') 1402 1403 write(message,'(3x,a,7x,e16.9,1x,a,i4,a,i4,a,a)')& 1404 & '+correction',polb(isppol),'(isppol=',isppol,', istep=1..',istep,')',ch10 1405 call wrtout(std_out,message,'COLL') 1406 1407 else 1408 1409 write(message,'(9x,a,7x,e16.9,1x,a,i4,a)')& 1410 & 'total',polb_mult(isppol,istep),'(isppol=',isppol,')'!,ch10 1411 call wrtout(std_out,message,'COLL') 1412 end if 1413 1414 polbtot = polbtot + polb(isppol) 1415 1416 end do ! isppol 1417 1418 ! Fold into interval [-1,1] 1419 polbtot = polbtot - 2_dp*nint(polbtot/2_dp) 1420 1421 ABI_FREE(det_string) 1422 ABI_FREE(polberry) 1423 1424 ! ========================================================================== 1425 1426 ! Compute the ionic Berry phase 1427 1428 call xred2xcart(natom,rprimd,xcart,xred) 1429 politot = zero 1430 write(message,'(a)')' Compute the ionic contributions' 1431 call wrtout(std_out,message,'COLL') 1432 1433 write(message,'(a,2x,a,2x,a,15x,a)')ch10,& 1434 & 'itom', 'itypat', 'polion' 1435 call wrtout(std_out,message,'COLL') 1436 1437 do iatom = 1, natom 1438 itypat = typat(iatom) 1439 1440 ! The ionic phase can be computed much easier 1441 polion = zion(itypat)*xred(idir,iatom) 1442 1443 ! Fold into interval (-1,1) 1444 polion = polion - 2_dp*nint(polion/2_dp) 1445 politot = politot + polion 1446 write(message,'(2x,i2,5x,i2,10x,e16.9)') iatom,itypat,polion 1447 call wrtout(std_out,message,'COLL') 1448 end do 1449 1450 ! Fold into interval [-1,1] again 1451 politot = politot - 2_dp*nint(politot/2_dp) 1452 pion(idir) = politot 1453 1454 write(message,'(9x,a,7x,es19.9)') 'total',politot 1455 call wrtout(std_out,message,'COLL') 1456 1457 1458 ! ========================================================================== 1459 1460 ! Compute the total polarization 1461 1462 poltot = politot + polbtot 1463 1464 if (berrystep==1)then 1465 write(message,'(a,a)')ch10,& 1466 & ' Summary of the results' 1467 call wrtout(std_out,message,'COLL') 1468 if (unit_out /= 0) then 1469 call wrtout(unit_out,message,'COLL') 1470 end if 1471 else 1472 write(message,'(a,a,i4)')ch10,& 1473 & ' Summary of the results for istep =',istep 1474 call wrtout(std_out,message,'COLL') 1475 if (unit_out /= 0) then 1476 call wrtout(unit_out,message,'COLL') 1477 end if 1478 end if 1479 1480 write(message,'(a,es19.9)')& 1481 & ' Electronic Berry phase ' ,polbtot 1482 call wrtout(std_out,message,'COLL') 1483 if (unit_out /= 0) then 1484 call wrtout(unit_out,message,'COLL') 1485 end if 1486 1487 write(message,'(a,es19.9)') & 1488 & ' Ionic phase ', politot 1489 call wrtout(std_out,message,'COLL') 1490 if (unit_out /= 0) then 1491 call wrtout(unit_out,message,'COLL') 1492 end if 1493 1494 write(message,'(a,es19.9)') & 1495 & ' Total phase ', poltot 1496 call wrtout(std_out,message,'COLL') 1497 if (unit_out /= 0) then 1498 call wrtout(unit_out,message,'COLL') 1499 end if 1500 1501 ! REC start 1502 if(abs(dtset%polcen(idir))>tol8)then 1503 poltot = poltot-dtset%polcen(idir) 1504 write(message,'(a,f15.10)') & 1505 & ' Translating Polarization by P0 for centrosymmetric cell: ',& 1506 & dtset%polcen(idir) 1507 call wrtout(std_out,message,'COLL') 1508 if (unit_out /= 0) then 1509 call wrtout(unit_out,message,'COLL') 1510 end if 1511 end if 1512 ! REC end 1513 1514 poltot = poltot - 2.0_dp*nint(poltot/2._dp) 1515 write(message,'(a,es19.9)') & 1516 & ' Remapping in [-1,1] ', poltot 1517 call wrtout(std_out,message,'COLL') 1518 if (unit_out /= 0) then 1519 call wrtout(unit_out,message,'COLL') 1520 end if 1521 1522 ! ! REC and HONG 1523 ! ===================================================================================== 1524 ! Polarization branch control (start) 1525 ! ------------------------------------------------------------------------------------- 1526 ! berrysav == 0, for non fixed D/d calculation, polarizaion is in [-1,1],done above 1527 ! for fixed D/d calculation, choose polarization to minimize internal 1528 ! energy, or minimize |red_efiled|. (red_dfield=red_efiled+red_ptot) 1529 ! (d=e+p, as (26) of Stengel, Suppl.) [[cite:Stengel2009]] 1530 ! This is default value. 1531 ! 1532 ! berrysav == 1, keep the polarization on the same branch, which saved in file POLSAVE 1533 ! ====================================================================================== 1534 1535 ! for fixed D/d calculation, choose polarization to minimize internal energy, or to minimize reduced electric field |red_efield| 1536 if(dtset%berrysav ==0 .and. (dtset%berryopt == 6 .or. dtset%berryopt == 7 .or. & 1537 & dtset%berryopt == 16 .or. dtset%berryopt == 17)) then 1538 1539 jump=-nint(dtset%red_dfield(idir) - poltot) ! red_efield = red_dfield - poltot 1540 1541 if(jump /= 0)then 1542 write(message,'(a,i1,a,es19.9,a,i2)') & 1543 & ' P(',idir,') Shifted polarization branch to minimize red_efield & 1544 & k from ',poltot, ' by ',jump 1545 call wrtout(std_out,message,'COLL') 1546 if (unit_out /= 0) then 1547 call wrtout(unit_out,message,'COLL') 1548 end if 1549 poltot=poltot-jump 1550 end if 1551 pol0(idir)=poltot 1552 1553 end if 1554 1555 1556 ! keep the polarization on the same branch. 1557 if (dtset%berrysav == 1) then 1558 1559 ! use saved polarization to keep on same branch 1560 inquire(file='POLSAVE',exist=lexist) 1561 if(lexist)then 1562 if(idir==1)then 1563 if(mpi_enreg%me==0)then 1564 if (open_file('POLSAVE',message,newunit=unt,status='OLD') /= 0) then 1565 ABI_ERROR(message) 1566 end if 1567 read(unt,*)pol0 1568 write(message,'(a,3f20.12)')'Reading old polarization:',pol0 1569 call wrtout(std_out,message,'COLL') 1570 if (unit_out /= 0) then 1571 call wrtout(unit_out,message,'COLL') 1572 end if 1573 close(unt) 1574 end if 1575 call xmpi_bcast(pol0,0,spaceComm,ierr) 1576 end if 1577 else 1578 pol0(idir)=poltot 1579 end if 1580 jump=nint(poltot-pol0(idir)) 1581 if(jump /= 0)then 1582 write(message,'(a,i1,a,es19.9,a,i2)') & 1583 & ' P(',idir,') jumped to new branch. Shifting bac& 1584 & k from ',poltot, ' by ',jump 1585 call wrtout(std_out,message,'COLL') 1586 if (unit_out /= 0) then 1587 call wrtout(unit_out,message,'COLL') 1588 end if 1589 poltot=poltot-jump 1590 end if 1591 1592 pol0(idir)=poltot 1593 1594 end if 1595 1596 ! ===================================================================================== 1597 ! Polarization branch control (end) 1598 ! ===================================================================================== 1599 1600 1601 ! Transform the phase into a polarization 1602 fac = 1._dp/(gmod*dtefield%nkstr(idir)) 1603 ! !REC fac = fac/ucvol 1604 ! !REC pol = fac*poltot 1605 red_ptot(idir)=poltot !!REC 1606 pol = fac*red_ptot(idir)/ucvol !!REC 1607 ptot(idir)=red_ptot(idir)/ucvol !!REC 1608 write(message,'(a,a,es19.9,a,a,a,es19.9,a,a)')ch10,& 1609 & ' Polarization ', pol,' (a.u. of charge)/bohr^2',ch10,& 1610 & ' Polarization ', pol*(e_Cb)/(Bohr_Ang*1d-10)**2,& 1611 & ' C/m^2',ch10 1612 call wrtout(std_out,message,'COLL') 1613 if (unit_out /= 0) then 1614 call wrtout(unit_out,message,'COLL') 1615 end if 1616 1617 1618 ABI_FREE(idxkstr_mult) 1619 1620 end do !istep 1621 1622 pel(idir) = polbtot 1623 1624 end if ! if calculate polarization polflag==1 1625 1626 end do ! Close loop over idir 1627 1628 !!REC start 1629 if (dtset%berrysav == 1) then 1630 if(mpi_enreg%me==0)then 1631 if (open_file('POLSAVE',message,newunit=unt,status='UNKNOWN') /= 0) then 1632 ABI_ERROR(message) 1633 end if 1634 write(unt,'(3F20.12)') pol0 1635 close(unt) 1636 end if 1637 first=.false. 1638 end if 1639 !!REC end 1640 1641 !------------------------------------------------- 1642 ! Compute polarization in cartesian coordinates 1643 !------------------------------------------------- 1644 if (all(dtset%rfdir(:) == 1)) then 1645 1646 if(usepaw.ne.1) then 1647 pelev=zero 1648 else 1649 call pawpolev(my_natom,natom,ntypat,pawrhoij,pawtab,pelev,& 1650 & comm_atom=mpi_enreg%comm_atom) 1651 ! note that in the PAW case, the pelev contribution is already 1652 ! implicitly included in the electronic polarization, from the 1653 ! discretized derivative operator. In the NCPP case no such 1654 ! terms exist anyway. Actually in the PAW formulation 1655 ! such terms are included to all orders, unlike in USPP where only 1656 ! zeroth and first-order terms are. In USPP the first-order term 1657 ! is pelev. Here we compute pelev separately only for reporting 1658 ! purposes in polcart, it is not added into pel or used in the the 1659 ! PAW finite field code in make_grad_berry.F90 1660 ! 13 June 2012 J Zwanziger 1661 end if 1662 call polcart(red_ptot,pel,pel_cart,pelev,pion,pion_cart,3,& 1663 & ptot_cart,rprimd,ucvol,unit_out,usepaw) 1664 1665 end if 1666 1667 !deallocate(pwind_k, dtm) 1668 ABI_FREE(pwnsfac_k) 1669 ABI_FREE(sflag_k) 1670 ABI_FREE(cg1_k) 1671 if (ddkflag == 1) then 1672 ABI_FREE(cg1) 1673 ABI_FREE(eig_dum) 1674 ABI_FREE(occ_dum) 1675 end if 1676 1677 if (usepaw == 1) then 1678 ABI_FREE(dimlmn) 1679 call pawcprj_free(cprj_k) 1680 call pawcprj_free(cprj_kb) 1681 call pawcprj_free(cprj_gat) 1682 ABI_FREE(cprj_k) 1683 ABI_FREE(cprj_kb) 1684 ABI_FREE(cprj_gat) 1685 if (dtset%kptopt /= 3) then 1686 call pawcprj_free(cprj_ikn) 1687 call pawcprj_free(cprj_fkn) 1688 ABI_FREE(cprj_ikn) 1689 ABI_FREE(cprj_fkn) 1690 end if 1691 if (calc_epaw3_force) then 1692 ABI_FREE(dsdr_sum) 1693 ABI_FREE(epawf3_str) 1694 end if 1695 if (calc_epaw3_stress) then 1696 ABI_FREE(dsds_sum) 1697 ABI_FREE(epaws3_str) 1698 end if 1699 1700 if (nproc>1) then 1701 call pawcprj_free(cprj_buf) 1702 ABI_FREE(cprj_buf) 1703 end if 1704 1705 end if 1706 1707 ABI_FREE(ikpt3) 1708 ABI_FREE(ikpt3i) 1709 ABI_FREE(sflag_k_mult) 1710 ABI_FREE(npw_k3) 1711 ABI_FREE(pwind_k_mult) 1712 ABI_FREE(itrs_mult) 1713 ABI_FREE(coef) 1714 ABI_FREE(polb_mult) 1715 ABI_FREE(dtm_mult) 1716 1717 !DEBUG 1718 !write(std_out,*)'berryphase_new exit' 1719 !END_DEBUG 1720 1721 end subroutine berryphase_new
ABINIT/init_e_field_vars [ Functions ]
NAME
init_e_field_vars
FUNCTION
Initialization of variables and data structures used in polarization calculations
COPYRIGHT
Copyright (C) 2004-2024 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
SOURCE
2713 subroutine init_e_field_vars(dtefield,dtset,gmet,gprimd,kg,& 2714 & mpi_enreg,npwarr,occ,pawang,pawrad,pawtab,psps,& 2715 & pwind,pwind_alloc,pwnsfac,rprimd,symrec,xred) 2716 2717 !Arguments ------------------------------------ 2718 !scalars 2719 integer,intent(out) :: pwind_alloc 2720 type(MPI_type),intent(inout) :: mpi_enreg 2721 type(dataset_type),intent(inout) :: dtset 2722 type(efield_type),intent(inout) :: dtefield !vz_i needs efield2 2723 type(pawang_type),intent(in) :: pawang 2724 type(pseudopotential_type),intent(in) :: psps 2725 !arrays 2726 integer,intent(in) :: kg(3,dtset%mpw*dtset%mkmem),npwarr(dtset%nkpt) 2727 integer,intent(in) :: symrec(3,3,dtset%nsym) 2728 integer,pointer :: pwind(:,:,:) 2729 real(dp),intent(in) :: gmet(3,3),gprimd(3,3),occ(dtset%mband*dtset%nkpt*dtset%nsppol) 2730 real(dp),intent(in) :: rprimd(3,3),xred(3,dtset%natom) 2731 real(dp),pointer :: pwnsfac(:,:) 2732 type(pawrad_type),intent(in) :: pawrad(dtset%ntypat*psps%usepaw) 2733 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*psps%usepaw) 2734 2735 !Local variables------------------------------- 2736 logical :: initfield 2737 !scalars 2738 2739 ! ************************************************************************* 2740 2741 initfield = .false. 2742 2743 !initialization 2744 dtefield%has_qijb = 0 2745 dtefield%has_epawf3 = 0 2746 dtefield%has_epaws3 = 0 2747 dtefield%has_expibi = 0 2748 dtefield%has_rij = 0 2749 dtefield%usecprj = 0 2750 dtefield%berryopt = 0 2751 2752 if ((dtset%berryopt < 0).or.(dtset%berryopt == 4) .or. (dtset%berryopt == 6) .or.(dtset%berryopt == 7) .or. & 2753 & (dtset%berryopt == 14) .or.(dtset%berryopt == 16) .or.(dtset%berryopt == 17)) then 2754 nullify(pwind,pwnsfac) 2755 call initberry(dtefield,dtset,gmet,gprimd,kg,& 2756 & dtset%mband,dtset%mkmem,mpi_enreg,dtset%mpw,& 2757 & dtset%natom,dtset%nkpt,npwarr,dtset%nsppol,& 2758 & dtset%nsym,dtset%ntypat,occ,pawang,pawrad,pawtab,& 2759 & psps,pwind,pwind_alloc,pwnsfac,rprimd,symrec,& 2760 & dtset%typat,psps%usepaw,xred) 2761 initfield = .true. 2762 end if 2763 2764 if (.not. initfield .and. dtset%orbmag == 0) then 2765 ! initorbmag.F90 also allocates pwind and pwnsfac 2766 pwind_alloc = 1 2767 ABI_MALLOC(pwind,(pwind_alloc,2,3)) 2768 ABI_MALLOC(pwnsfac,(2,pwind_alloc)) 2769 pwind(:,:,:)=0 2770 pwnsfac(:,:)=zero 2771 end if 2772 2773 end subroutine init_e_field_vars
ABINIT/initberry [ 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-2024 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 = information 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)
SOURCE
2849 subroutine initberry(dtefield,dtset,gmet,gprimd,kg,mband,& 2850 & mkmem,mpi_enreg,mpw,natom,nkpt,npwarr,nsppol,& 2851 & nsym,ntypat,occ,pawang,pawrad,pawtab,psps,& 2852 & pwind,pwind_alloc,pwnsfac,& 2853 & rprimd,symrec,typat,usepaw,xred) 2854 2855 !Arguments ------------------------------------ 2856 !scalars 2857 integer,intent(in) :: mband,mkmem,mpw,natom,nkpt,nsppol,nsym,ntypat,usepaw 2858 integer,intent(out) :: pwind_alloc 2859 type(MPI_type),intent(inout) :: mpi_enreg 2860 type(dataset_type),intent(inout) :: dtset 2861 type(efield_type),intent(inout) :: dtefield !vz_i 2862 type(pawang_type),intent(in) :: pawang 2863 type(pseudopotential_type),intent(in) :: psps 2864 !arrays 2865 integer,intent(in) :: kg(3,mpw*mkmem),npwarr(nkpt) 2866 integer,intent(in) :: symrec(3,3,nsym),typat(natom) 2867 integer,pointer :: pwind(:,:,:) 2868 real(dp),intent(in) :: gmet(3,3),gprimd(3,3),occ(mband*nkpt*nsppol) 2869 real(dp),intent(in) :: rprimd(3,3),xred(3,natom) 2870 real(dp),pointer :: pwnsfac(:,:) 2871 type(pawrad_type),intent(in) :: pawrad(ntypat) 2872 type(pawtab_type),intent(in) :: pawtab(ntypat) 2873 2874 !Local variables------------------------------- 2875 !scalars 2876 integer :: exchn2n3d,flag,flag_kpt,fnkpt_computed,iband,icg,icprj 2877 integer :: idir,idum,idum1,ierr,ifor,ikg,ikg1,ikpt,ikpt1,ikpt1f 2878 integer :: ikpt1i,ikpt2,ikpt_loc,ikptf,ikpti,ikstr,index,ineigh,ipw,ipwnsfac 2879 integer :: isppol,istr,istwf_k,isym,isym1,itrs,itypat,iunmark,jpw,klmn,lmax,lmn2_size_max 2880 integer :: me,me_g0,mkmem_,my_nspinor,nband_k,mband_occ_k,ncpgr,nkstr,nproc,npw_k,npw_k1,spaceComm 2881 integer :: option, brav, mkpt, nkptlatt 2882 integer :: jstr,ii,jj,isign 2883 integer :: dk_flag, coord1, coord2 2884 integer :: mult 2885 real(dp) :: c1,ecut_eff,eg,eg_ev,rdum,diffk1,diffk2,diffk3 2886 real(dp) :: dist_, max_dist, last_dist, dist,kpt_shifted1,kpt_shifted2,kpt_shifted3 2887 real(dp) :: gprimdlc(3,3),rmetllc(3,3),gmetllc(3,3),ucvol_local 2888 ! gprimd(3,3) = inverse of rprimd 2889 ! rmetlcl(3,3)=real-space metric (same as rmet in metric.F90) 2890 ! gmetlcl(3,3)= same as gmet in metric.F90 2891 ! ucvol = volume of the unit cell in Bohr**3 2892 character(len=500) :: message 2893 logical :: calc_epaw3_force,calc_epaw3_stress,fieldflag 2894 !arrays 2895 integer :: dg(3),iadum(3),iadum1(3),neigh(6) 2896 integer,allocatable :: dimlmn(:),kg1_k(:,:),kpt_mark(:),nattyp_dum(:) 2897 real(dp) :: diffk(3),dk(3),dum33(3,3),eg_dir(3) 2898 real(dp) :: kpt1(3) 2899 real(dp) :: delta_str3(2), dstr(2),dk_str(2,2,3) 2900 real(dp) :: tsec(2) 2901 real(dp),allocatable :: calc_expibi(:,:),calc_qijb(:,:,:),spkpt(:,:) 2902 2903 ! ************************************************************************* 2904 2905 DBG_ENTER("COLL") 2906 2907 call timab(1001,1,tsec) 2908 call timab(1002,1,tsec) 2909 2910 spaceComm=mpi_enreg%comm_cell 2911 nproc=xmpi_comm_size(spaceComm) 2912 me=xmpi_comm_rank(spaceComm) 2913 2914 !save the current value of berryopt 2915 dtefield%berryopt = dtset%berryopt 2916 !save the current value of nspinor 2917 dtefield%nspinor = dtset%nspinor 2918 2919 !---------------------------------------------------------------------------- 2920 !-------------------- Obtain k-point grid in the full BZ -------------------- 2921 !---------------------------------------------------------------------------- 2922 2923 if(dtset%kptopt==1 .or. dtset%kptopt==2 .or. dtset%kptopt==4)then 2924 ! Compute the number of k points in the G-space unit cell 2925 nkptlatt=dtset%kptrlatt(1,1)*dtset%kptrlatt(2,2)*dtset%kptrlatt(3,3) & 2926 & +dtset%kptrlatt(1,2)*dtset%kptrlatt(2,3)*dtset%kptrlatt(3,1) & 2927 & +dtset%kptrlatt(1,3)*dtset%kptrlatt(2,1)*dtset%kptrlatt(3,2) & 2928 & -dtset%kptrlatt(1,2)*dtset%kptrlatt(2,1)*dtset%kptrlatt(3,3) & 2929 & -dtset%kptrlatt(1,3)*dtset%kptrlatt(2,2)*dtset%kptrlatt(3,1) & 2930 & -dtset%kptrlatt(1,1)*dtset%kptrlatt(2,3)*dtset%kptrlatt(3,2) 2931 2932 ! Call smpbz to obtain the list of k-point in the full BZ - without symmetry reduction 2933 option = 0 2934 brav = 1 2935 mkpt=nkptlatt*dtset%nshiftk 2936 ABI_MALLOC(spkpt,(3,mkpt)) 2937 call smpbz(1,ab_out,dtset%kptrlatt,mkpt,fnkpt_computed,dtset%nshiftk,option,dtset%shiftk,spkpt) 2938 dtefield%fnkpt = fnkpt_computed 2939 ABI_MALLOC(dtefield%fkptns,(3,dtefield%fnkpt)) 2940 dtefield%fkptns(:,:)=spkpt(:,1:dtefield%fnkpt) 2941 ABI_FREE(spkpt) 2942 else if(dtset%kptopt==3.or.dtset%kptopt==0)then 2943 dtefield%fnkpt=nkpt 2944 ABI_MALLOC(dtefield%fkptns,(3,dtefield%fnkpt)) 2945 dtefield%fkptns(1:3,1:dtefield%fnkpt)=dtset%kpt(1:3,1:dtefield%fnkpt) 2946 if(dtset%kptopt==0)then 2947 write(message,'(10a)') ch10,& 2948 & ' initberry : WARNING -',ch10,& 2949 & ' you have defined manually the k-point grid with kptopt = 0',ch10,& 2950 & ' the berry phase calculation works only with a regular k-points grid,',ch10,& 2951 & ' abinit doesn''t check if your grid is regular...' 2952 call wrtout(std_out,message,'PERS') 2953 end if 2954 end if 2955 2956 !call listkk to get mapping from FBZ to IBZ 2957 rdum=1.0d-5 ! cutoff distance to decide when two k points match 2958 ABI_MALLOC(dtefield%indkk_f2ibz,(dtefield%fnkpt,6)) 2959 2960 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 2961 2962 !ji: The following may need modification in the future 2963 !**** no spin-polarization doubling ; allow use of time reversal symmetry **** 2964 2965 !Here is original call 2966 ! 2967 !call listkk(rdum,gmet,dtefield%indkk_f2ibz,dtset%kptns,dtefield%fkptns,nkpt,& 2968 !& dtefield%fnkpt,dtset%nsym,1,dtset%symafm,dtset%symrel,1, spaceComm) 2969 2970 call timab(1002,2,tsec) 2971 call timab(1003,1,tsec) 2972 2973 call listkk(rdum,gmet,dtefield%indkk_f2ibz,dtset%kptns,dtefield%fkptns,nkpt,& 2974 & dtefield%fnkpt,dtset%nsym,1,dtset%symafm,symrec,1, spaceComm, use_symrec=.True.) 2975 2976 call timab(1003,2,tsec) 2977 call timab(1004,1,tsec) 2978 2979 !Construct i2fbz and f2ibz 2980 ABI_MALLOC(dtefield%i2fbz,(nkpt)) 2981 idum=0 2982 do ikpt=1,dtefield%fnkpt 2983 if (dtefield%indkk_f2ibz(ikpt,2)==1 .and. & 2984 & dtefield%indkk_f2ibz(ikpt,6) == 0 .and. & 2985 & maxval(abs(dtefield%indkk_f2ibz(ikpt,3:5))) == 0 ) then 2986 dtefield%i2fbz(dtefield%indkk_f2ibz(ikpt,1))=ikpt 2987 idum=idum+1 2988 end if 2989 end do 2990 if (idum/=nkpt)then 2991 ABI_ERROR('Found wrong number of k-points in IBZ') 2992 end if 2993 2994 !set flags for fields, forces, stresses 2995 fieldflag = ( (dtset%berryopt== 4) .or. (dtset%berryopt== 6) .or. (dtset%berryopt== 7) & 2996 & .or. (dtset%berryopt==14) .or. (dtset%berryopt==16) .or. (dtset%berryopt==17) ) 2997 ! following two flags activates computation of projector gradient contributions to force and 2998 ! stress in finite field PAW calculations 2999 calc_epaw3_force = (fieldflag .and. usepaw == 1 .and. dtset%optforces /= 0) 3000 calc_epaw3_stress = (fieldflag .and. usepaw == 1 .and. dtset%optstress /= 0) 3001 3002 3003 3004 !---------------------------------------------------------------------------- 3005 !------------- Allocate PAW space if necessary ------------------------------ 3006 !---------------------------------------------------------------------------- 3007 3008 if (usepaw == 1) then 3009 3010 dtefield%usepaw = usepaw 3011 dtefield%natom = natom 3012 dtefield%my_natom = mpi_enreg%my_natom 3013 3014 ABI_MALLOC(dtefield%lmn_size,(ntypat)) 3015 ABI_MALLOC(dtefield%lmn2_size,(ntypat)) 3016 do itypat = 1, ntypat 3017 dtefield%lmn_size(itypat) = pawtab(itypat)%lmn_size 3018 dtefield%lmn2_size(itypat) = pawtab(itypat)%lmn2_size 3019 end do 3020 3021 lmn2_size_max = psps%lmnmax*(psps%lmnmax+1)/2 3022 dtefield%lmn2max = lmn2_size_max 3023 3024 ! expibi and qijb_kk are NOT parallelized over atoms 3025 ! this may change in the future (JZwanziger 18 March 2014) 3026 ABI_MALLOC(dtefield%qijb_kk,(2,lmn2_size_max,dtefield%natom,3)) 3027 ABI_MALLOC(dtefield%expibi,(2,dtefield%natom,3)) 3028 dtefield%has_expibi = 1 3029 dtefield%has_qijb = 1 3030 3031 if ( fieldflag .and. dtefield%has_rij==0) then 3032 lmn2_size_max = psps%lmnmax*(psps%lmnmax+1)/2 3033 ABI_MALLOC(dtefield%rij,(lmn2_size_max,ntypat,3)) 3034 dtefield%has_rij = 1 3035 end if 3036 3037 ! additional F3-type force term for finite electric field with PAW. Same term 3038 ! might also apply for other displacement-type field calculations, but not sure yet 3039 ! JZwanziger 4 April 2014 3040 if ( calc_epaw3_force ) then 3041 ABI_MALLOC(dtefield%epawf3,(dtefield%natom,3,3)) 3042 dtefield%has_epawf3 = 1 3043 end if 3044 if ( calc_epaw3_stress ) then 3045 ABI_MALLOC(dtefield%epaws3,(dtefield%natom,3,6)) 3046 dtefield%has_epaws3 = 1 3047 end if 3048 3049 ncpgr = 0 3050 if ( fieldflag .and. dtefield%usecprj == 0) then 3051 ABI_MALLOC(dimlmn,(natom)) 3052 call pawcprj_getdim(dimlmn,natom,nattyp_dum,ntypat,typat,pawtab,'R') 3053 ! allocate space for cprj at kpts in BZ (IBZ or FBZ) 3054 ABI_MALLOC(dtefield%cprj,(natom, mband*dtset%nspinor*dtset%nkpt*nsppol)) 3055 ! write(std_out,*) "initberry alloc of cprj ", shape(dtefield%cprj) 3056 if (calc_epaw3_force .and. .not. calc_epaw3_stress) ncpgr = 3 3057 if (.not. calc_epaw3_force .and. calc_epaw3_stress) ncpgr = 6 3058 if (calc_epaw3_force .and. calc_epaw3_stress) ncpgr = 9 3059 call pawcprj_alloc(dtefield%cprj,ncpgr,dimlmn) 3060 dtefield%usecprj = 1 3061 ABI_FREE(dimlmn) 3062 end if 3063 3064 ABI_MALLOC(dtefield%cprjindex,(nkpt,nsppol)) 3065 dtefield%cprjindex(:,:) = 0 3066 3067 if (dtset%kptopt /= 3) then 3068 ABI_MALLOC(dtefield%atom_indsym,(4,nsym,natom)) 3069 call symatm(dtefield%atom_indsym,natom,nsym,symrec,dtset%tnons,tol8,typat,xred) 3070 lmax = psps%mpsang - 1 3071 ABI_MALLOC(dtefield%zarot,(2*lmax+1,2*lmax+1,lmax+1,nsym)) 3072 call setsym_ylm(gprimd,lmax,nsym,1,rprimd,symrec,dtefield%zarot) 3073 dtefield%nsym = nsym 3074 dtefield%lmax = lmax 3075 dtefield%lmnmax = psps%lmnmax 3076 end if 3077 3078 end if 3079 3080 !------------------------------------------------------------------------------ 3081 !------------------- Compute variables related to MPI // ---------------------- 3082 !------------------------------------------------------------------------------ 3083 3084 3085 3086 3087 if (nproc==1) then 3088 dtefield%fmkmem = dtefield%fnkpt 3089 dtefield%fmkmem_max = dtefield%fnkpt 3090 dtefield%mkmem_max = nkpt 3091 else 3092 dtefield%fmkmem = 0 3093 do ikpt = 1, dtefield%fnkpt 3094 ikpti = dtefield%indkk_f2ibz(ikpt,1) 3095 nband_k = dtset%nband(ikpti) 3096 if (.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpti,1,nband_k,-1,me))) & 3097 & dtefield%fmkmem = dtefield%fmkmem + 1 3098 end do 3099 ! Maximum value of mkmem and fmkmem 3100 call xmpi_max(dtefield%fmkmem,dtefield%fmkmem_max,spaceComm,ierr) 3101 ! I have to use the dummy variable mkmem_ because 3102 ! mkmem is declared as intent(in) while the first 3103 ! argument of xmpi_max must be intent(inout) 3104 mkmem_ = mkmem 3105 call xmpi_max(mkmem_,dtefield%mkmem_max,spaceComm,ierr) 3106 end if 3107 3108 ABI_MALLOC(mpi_enreg%kpt_loc2fbz_sp,(0:nproc-1,1:dtefield%fmkmem_max*nsppol, 1:2)) 3109 ABI_MALLOC(mpi_enreg%kpt_loc2ibz_sp,(0:nproc-1,1:dtefield%mkmem_max*nsppol, 1:2)) 3110 ABI_MALLOC(mpi_enreg%kptdstrb,(nproc,6,dtefield%fmkmem_max*nsppol*2)) 3111 ABI_MALLOC(mpi_enreg%mkmem,(0:nproc-1)) 3112 mpi_enreg%kpt_loc2fbz_sp(:,:,:) = 0 3113 mpi_enreg%kpt_loc2ibz_sp(:,:,:) = 0 3114 mpi_enreg%kptdstrb(:,:,:) = 0 3115 mpi_enreg%mkmem(:) = 0 3116 3117 if (fieldflag) then 3118 ABI_MALLOC(dtefield%cgqindex,(3,6,nkpt*nsppol)) 3119 ABI_MALLOC(dtefield%nneigh,(nkpt)) 3120 dtefield%cgqindex(:,:,:) = 0 ; dtefield%nneigh(:) = 0 3121 end if 3122 3123 pwind_alloc = mpw*dtefield%fmkmem_max 3124 ABI_MALLOC(pwind,(pwind_alloc,2,3)) 3125 ABI_MALLOC(pwnsfac,(2,pwind_alloc)) 3126 3127 !------------------------------------------------------------------------------ 3128 !---------------------- Compute efield_type variables ------------------------- 3129 !------------------------------------------------------------------------------ 3130 3131 !Initialization of efield_type variables 3132 mult=dtset%useria+1 3133 dtefield%efield_dot(:) = zero 3134 dtefield%dkvecs(:,:) = zero 3135 dtefield%maxnstr = 0 ; dtefield%maxnkstr = 0 3136 dtefield%nstr(:) = 0 ; dtefield%nkstr(:) = 0 3137 ABI_MALLOC(dtefield%ikpt_dk,(dtefield%fnkpt,2,3)) 3138 ABI_MALLOC(dtefield%cgindex,(nkpt,nsppol)) 3139 ABI_MALLOC(dtefield%kgindex,(nkpt)) 3140 ABI_MALLOC(dtefield%fkgindex,(dtefield%fnkpt)) 3141 dtefield%ikpt_dk(:,:,:) = 0 3142 dtefield%cgindex(:,:) = 0 3143 dtefield%mband_occ = 0 3144 ABI_MALLOC(dtefield%nband_occ,(nsppol)) 3145 dtefield%kgindex(:) = 0 3146 dtefield%fkgindex(:) = 0 3147 3148 if (fieldflag) then 3149 dtset%rfdir(1:3) = 1 3150 end if 3151 3152 3153 !Compute spin degeneracy 3154 if (nsppol == 1 .and. dtset%nspinor == 1) then 3155 dtefield%sdeg = two 3156 else if (nsppol == 2 .or. my_nspinor == 2) then 3157 dtefield%sdeg = one 3158 end if 3159 3160 !Compute the number of occupied bands and check that 3161 !it is the same for each k-point 3162 3163 index = 0 3164 do isppol = 1, nsppol 3165 dtefield%nband_occ(isppol) = 0 3166 do ikpt = 1, nkpt 3167 3168 mband_occ_k = 0 3169 nband_k = dtset%nband(ikpt + (isppol - 1)*nkpt) 3170 3171 do iband = 1, nband_k 3172 index = index + 1 3173 if (abs(occ(index) - dtefield%sdeg) < tol8) mband_occ_k = mband_occ_k + 1 3174 end do 3175 3176 if (fieldflag) then 3177 if (nband_k /= mband_occ_k) then 3178 write(message,'(a,a,a)')& 3179 & ' In a finite electric field, nband must be equal ',ch10,& 3180 & ' to the number of valence bands.' 3181 ABI_ERROR(message) 3182 end if 3183 end if 3184 3185 if (ikpt > 1) then 3186 if (dtefield%nband_occ(isppol) /= mband_occ_k) then 3187 message = "The number of valence bands is not the same for every k-point of present spin channel" 3188 ABI_ERROR(message) 3189 end if 3190 else 3191 dtefield%mband_occ = max(dtefield%mband_occ, mband_occ_k) 3192 dtefield%nband_occ(isppol) = mband_occ_k 3193 end if 3194 3195 end do ! close loop over ikpt 3196 end do ! close loop over isppol 3197 3198 if (fieldflag) then 3199 ABI_MALLOC(dtefield%smat,(2,dtefield%mband_occ,dtefield%mband_occ,nkpt*nsppol,2,3)) 3200 3201 dtefield%smat(:,:,:,:,:,:) = zero 3202 end if 3203 3204 ABI_MALLOC(dtefield%sflag,(dtefield%mband_occ,nkpt*nsppol,2,3)) 3205 dtefield%sflag(:,:,:,:) = 0 3206 3207 !Compute the location of each wavefunction 3208 3209 icg = 0 3210 icprj = 0 3211 !ikg = 0 3212 do isppol = 1, nsppol 3213 do ikpt = 1, nkpt 3214 3215 nband_k = dtset%nband(ikpt + (isppol-1)*nkpt) 3216 3217 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) cycle 3218 3219 dtefield%cgindex(ikpt,isppol) = icg 3220 npw_k = npwarr(ikpt) 3221 icg = icg + npw_k*dtefield%nspinor*nband_k 3222 3223 if (usepaw == 1) then 3224 dtefield%cprjindex(ikpt,isppol) = icprj 3225 icprj = icprj + dtefield%nspinor*nband_k 3226 end if 3227 3228 end do 3229 end do 3230 3231 ikg = 0 3232 do ikpt = 1, nkpt 3233 if ((proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,1,me)).and.& 3234 & (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,nsppol,me))) cycle 3235 3236 npw_k = npwarr(ikpt) 3237 dtefield%kgindex(ikpt) = ikg 3238 ikg = ikg + npw_k 3239 end do 3240 3241 !Need to use dtset%red_efieldbar in the whole code 3242 !Compute the reciprocal lattice coordinates of the electric field 3243 if (fieldflag) then 3244 3245 call metric(gmetllc,gprimdlc,-1,rmetllc,rprimd,ucvol_local) 3246 3247 if (dtset%berryopt == 4 .or. dtset%berryopt == 6 .or. dtset%berryopt == 7) then 3248 3249 do ii=1,3 3250 dtset%red_efieldbar(ii) = dot_product(dtset%efield(:),rprimd(:,ii)) 3251 dtefield%efield_dot(ii) = dtset%red_efieldbar(ii) 3252 end do 3253 3254 ! dtefield%efield_dot(1) = dot_product(dtset%efield(:),rprimd(:,1)) 3255 ! dtefield%efield_dot(2) = dot_product(dtset%efield(:),rprimd(:,2)) 3256 ! dtefield%efield_dot(3) = dot_product(dtset%efield(:),rprimd(:,3)) 3257 3258 write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,& 3259 & ' initberry: Reduced electric field (ebar)',ch10,& 3260 & ' red_efieldbar(1:3) = ',dtset%red_efieldbar(1:3),ch10 3261 call wrtout(std_out,message,'COLL') 3262 3263 end if 3264 3265 if (dtset%berryopt == 6 .or. dtset%berryopt ==7 ) then 3266 3267 do ii=1,3 3268 dtset%red_dfield(ii)= (dot_product(dtset%dfield(:),gprimdlc(:,ii)))*ucvol_local/(4.d0*pi) 3269 end do 3270 3271 write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,& 3272 & ' initberry: Reduced electric displacement field',ch10,& 3273 & ' red_dfield(1:3) = ',dtset%red_dfield(1:3),ch10 3274 call wrtout(std_out,message,'COLL') 3275 3276 end if 3277 3278 3279 if ( dtset%berryopt == 14 ) then 3280 ! transfer to unreduced electric field. 3281 do idir=1,3 3282 dtset%efield(idir)= dot_product(dtset%red_efieldbar(:),gprimdlc(:,idir)) 3283 dtefield%efield_dot(idir) = dtset%red_efieldbar(idir) 3284 ! dtefield%efield2(idir)=dtset%red_efieldbar(idir) 3285 end do 3286 3287 ! dtefield%efield_dot(1) = dtset%red_efieldbar(1) 3288 ! dtefield%efield_dot(2) = dtset%red_efieldbar(2) 3289 ! dtefield%efield_dot(3) = dtset%red_efieldbar(3) 3290 3291 write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,& 3292 & ' initberry: Unreduced electric field (a.u.)',ch10,& 3293 & ' efield(1:3) = ',dtset%efield(1:3),ch10 3294 call wrtout(std_out,message,'COLL') 3295 3296 write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,& 3297 & ' initberry: Reduced electric field (ebar)',ch10,& 3298 & ' red_efieldbar(1:3) = ',dtset%red_efieldbar(1:3),ch10 3299 call wrtout(std_out,message,'COLL') 3300 3301 end if 3302 3303 3304 if ( dtset%berryopt == 16 ) then 3305 3306 ! to calculate D 3307 do ii=1,3 3308 dtset%dfield(ii) =(4*pi/ucvol_local)*dot_product(dtset%red_dfield(:),rprimd(:,ii)) 3309 end do 3310 3311 do idir=1,3 3312 dtset%efield(idir)= (4*pi/ucvol_local)*dot_product(dtset%red_efield(:),rprimd(:,idir)) 3313 end do 3314 3315 do idir=1,3 3316 dtset%red_efieldbar(idir)= (4*pi/ucvol_local)*dot_product(dtset%red_efield(:),rmetllc(:,idir)) 3317 dtefield%efield_dot(idir) = dtset%red_efieldbar(idir) 3318 end do 3319 3320 ! dtefield%efield_dot(1) = dtset%red_efieldbar(1) 3321 ! dtefield%efield_dot(2) = dtset%red_efieldbar(2) 3322 ! dtefield%efield_dot(3) = dtset%red_efieldbar(3) 3323 3324 3325 write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,& 3326 & ' initberry: Unreduced electric displacement field (a.u.)',ch10,& 3327 & ' dfield(1:3) = ',dtset%dfield(1:3),ch10 3328 call wrtout(std_out,message,'COLL') 3329 3330 write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,& 3331 & ' initberry: Unreduced electric field (a.u.)',ch10,& 3332 & ' efield(1:3) = ',dtset%efield(1:3),ch10 3333 call wrtout(std_out,message,'COLL') 3334 3335 write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,& 3336 & ' initberry: Reduced electric field (ebar)',ch10,& 3337 & ' red_efieldbar(1:3) = ',dtset%red_efieldbar(1:3),ch10 3338 call wrtout(std_out,message,'COLL') 3339 3340 end if 3341 3342 if ( dtset%berryopt ==17) then 3343 3344 ! to calculate D 3345 3346 do idir=1,3 3347 dtset%efield(idir)= dot_product(dtset%red_efieldbar(:),gprimdlc(:,idir)) ! from ebar 3348 dtset%dfield(idir) =(4*pi/ucvol_local)*dot_product(dtset%red_dfield(:),rprimd(:,idir)) 3349 ! dtset%red_efield(idir) = (ucvol_local/(4*pi))*dot_product(dtset%red_efieldbar(:),gmetllc(:,idir)) 3350 dtefield%efield_dot(idir) = dtset%red_efieldbar(idir) 3351 end do 3352 3353 write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,& 3354 & ' initberry: Reduced electric field (ebar)',ch10,& 3355 & ' red_efieldbar(1:3) = ',dtset%red_efieldbar(1:3),ch10 3356 call wrtout(std_out,message,'COLL') 3357 3358 3359 write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,& 3360 & ' initberry: Unreduced electric field (a.u.)',ch10,& 3361 & ' efield(1:3) = ',dtset%efield(1:3),ch10 3362 call wrtout(std_out,message,'COLL') 3363 3364 write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,& 3365 & ' initberry: Reduced electric displacement field (a.u.)',ch10,& 3366 & ' red_dfield(1:3) = ',dtset%red_dfield(1:3),ch10 3367 call wrtout(std_out,message,'COLL') 3368 3369 write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,& 3370 & ' initberry: Unreduced electric displacement field (a.u.)',ch10,& 3371 & ' dfield(1:3) = ',dtset%dfield(1:3),ch10 3372 call wrtout(std_out,message,'COLL') 3373 3374 3375 end if 3376 3377 3378 3379 end if 3380 3381 call timab(1004,2,tsec) 3382 3383 !------------------------------------------------------------------------------ 3384 !---------------------- Compute dk -------------------------------------------- 3385 !------------------------------------------------------------------------------ 3386 3387 call timab(1005,1,tsec) 3388 3389 do idir = 1, 3 3390 3391 if (dtset%rfdir(idir) == 1) then 3392 3393 ! Compute dk(:), the vector between a k-point and its nearest 3394 ! neighbour along the direction idir 3395 3396 dk(:) = zero 3397 dk(idir) = 1._dp ! 1 mean there is no other k-point un the direction idir 3398 do ikpt = 2, dtefield%fnkpt 3399 diffk(:) = abs(dtefield%fkptns(:,ikpt) - dtefield%fkptns(:,1)) 3400 if ((diffk(1) < dk(1)+tol8).and.(diffk(2) < dk(2)+tol8).and.& 3401 & (diffk(3) < dk(3)+tol8)) dk(:) = diffk(:) 3402 end do 3403 dtefield%dkvecs(:,idir) = dk(:) 3404 ! DEBUG 3405 ! write(std_out,*)' initberry : idir, dk', idir, dk 3406 ! ENDDEBUG 3407 3408 ! For each k point, find k_prim such that k_prim= k + dk mod(G) 3409 ! where G is a vector of the reciprocal lattice 3410 3411 do ikpt = 1, dtefield%fnkpt 3412 3413 ! First k+dk, then k-dk 3414 do isign=-1,1,2 3415 kpt_shifted1=dtefield%fkptns(1,ikpt)- isign*dk(1) 3416 kpt_shifted2=dtefield%fkptns(2,ikpt)- isign*dk(2) 3417 kpt_shifted3=dtefield%fkptns(3,ikpt)- isign*dk(3) 3418 ! Note that this is still a order fnkpt**2 algorithm. 3419 ! It is possible to implement a order fnkpt algorithm, see listkk.F90. 3420 do ikpt1 = 1, dtefield%fnkpt 3421 diffk1=dtefield%fkptns(1,ikpt1) - kpt_shifted1 3422 if(abs(diffk1-nint(diffk1))>tol8)cycle 3423 diffk2=dtefield%fkptns(2,ikpt1) - kpt_shifted2 3424 if(abs(diffk2-nint(diffk2))>tol8)cycle 3425 diffk3=dtefield%fkptns(3,ikpt1) - kpt_shifted3 3426 if(abs(diffk3-nint(diffk3))>tol8)cycle 3427 dtefield%ikpt_dk(ikpt,(isign+3)/2,idir) = ikpt1 3428 exit 3429 end do ! ikpt1 3430 end do ! isign 3431 3432 ! OLD CODING 3433 ! First: k + dk 3434 ! do ikpt1 = 1, dtefield%fnkpt 3435 ! diffk(:) = abs(dtefield%fkptns(:,ikpt1) - & 3436 ! & dtefield%fkptns(:,ikpt) - dk(:)) 3437 ! if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then 3438 ! dtefield%ikpt_dk(ikpt,1,idir) = ikpt1 3439 ! exit 3440 ! end if 3441 ! end do 3442 3443 ! Second: k - dk 3444 ! do ikpt1 = 1, dtefield%fnkpt 3445 ! diffk(:) = abs(dtefield%fkptns(:,ikpt1) - & 3446 ! & dtefield%fkptns(:,ikpt) + dk(:)) 3447 ! if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then 3448 ! dtefield%ikpt_dk(ikpt,2,idir) = ikpt1 3449 ! exit 3450 ! end if 3451 ! end do 3452 3453 end do ! ikpt 3454 3455 ! Find the string length, starting from k point 1 3456 ! (all strings must have the same number of points) 3457 3458 nkstr = 1 3459 ikpt1 = 1 3460 do ikpt = 1, dtefield%fnkpt 3461 ikpt1 = dtefield%ikpt_dk(ikpt1,1,idir) 3462 if (ikpt1 == 1) exit 3463 nkstr = nkstr + 1 3464 end do 3465 3466 ! Check that the string length is a divisor of nkpt 3467 if(mod(dtefield%fnkpt,nkstr) /= 0) then 3468 write(message,'(a,i5,a,i7)')& 3469 & ' The string length = ',nkstr,& 3470 & ', is not a divisor of fnkpt =',dtefield%fnkpt 3471 ABI_BUG(message) 3472 end if 3473 3474 dtefield%nkstr(idir) = nkstr 3475 dtefield%nstr(idir) = dtefield%fnkpt/nkstr 3476 3477 end if ! dtset%rfdir(idir) == 1 3478 3479 write(message,'(a,i1,a,i3,a,i6)')& 3480 & ' initberry: for direction ',idir,', nkstr = ',dtefield%nkstr(idir),& 3481 & ', nstr = ',dtefield%nstr(idir) 3482 call wrtout(std_out,message,'COLL') 3483 call wrtout(ab_out,message,'COLL') 3484 3485 end do ! close loop over idir 3486 3487 call timab(1005,2,tsec) 3488 call timab(1006,1,tsec) 3489 3490 dtefield%maxnstr = maxval(dtefield%nstr(:)) 3491 dtefield%maxnkstr = maxval(dtefield%nkstr(:)) 3492 ABI_MALLOC(dtefield%idxkstr,(dtefield%maxnkstr,dtefield%maxnstr,3)) 3493 dtefield%idxkstr(:,:,:) = 0 3494 3495 !for the geometry of the string space : 3496 ABI_MALLOC(dtefield%coord_str,(2,dtefield%maxnstr,3)) 3497 ABI_MALLOC(dtefield%str_neigh,(-2:2,dtefield%maxnstr,3)) 3498 ABI_MALLOC(dtefield%strg_neigh,(-2:2,dtefield%maxnstr,2,3)) 3499 dtefield%coord_str(:,:,:) = 0.d0 3500 dtefield%str_neigh(:,:,:)=0 3501 dtefield%strg_neigh(:,:,:,:)=0 3502 dtefield%gmet_str(:,:,:)=0.d0 3503 3504 !------------------------------------------------------------------------------ 3505 !---------------------- Build the strings ------------------------------------- 3506 !------------------------------------------------------------------------------ 3507 3508 ABI_MALLOC(kpt_mark,(dtefield%fnkpt)) 3509 do idir = 1, 3 3510 3511 if (dtset%rfdir(idir) == 1) then 3512 3513 iunmark = 1 3514 kpt_mark(:) = 0 3515 do istr = 1, dtefield%nstr(idir) 3516 3517 do while(kpt_mark(iunmark) /= 0) 3518 iunmark = iunmark + 1 3519 end do 3520 dtefield%idxkstr(1,istr,idir) = iunmark 3521 kpt_mark(iunmark) = 1 3522 do ikstr = 2, dtefield%nkstr(idir) 3523 ikpt1 = dtefield%idxkstr(ikstr-1,istr,idir) 3524 ikpt2 = dtefield%ikpt_dk(ikpt1,1,idir) 3525 dtefield%idxkstr(ikstr,istr,idir) = ikpt2 3526 kpt_mark(ikpt2) = 1 3527 end do 3528 3529 end do ! istr 3530 3531 ! compute distance between strings 3532 ! compute the metric matrix of the strings space in the direction idir 3533 do ii = 1,3 3534 do jj = 1,3 3535 if (ii<idir.and.jj<idir) dtefield%gmet_str(ii ,jj ,idir) = & 3536 & gmet(ii,jj) - gmet(ii,idir)*gmet(jj,idir)/gmet(idir,idir) 3537 if (ii<idir.and.jj>idir) dtefield%gmet_str(ii ,jj-1,idir) = & 3538 & gmet(ii,jj) - gmet(ii,idir)*gmet(jj,idir)/gmet(idir,idir) 3539 if (ii>idir.and.jj<idir) dtefield%gmet_str(ii-1,jj ,idir) = & 3540 & gmet(ii,jj) - gmet(ii,idir)*gmet(jj,idir)/gmet(idir,idir) 3541 if (ii>idir.and.jj>idir) dtefield%gmet_str(ii-1,jj-1,idir) = & 3542 & gmet(ii,jj) - gmet(ii,idir)*gmet(jj,idir)/gmet(idir,idir) 3543 end do 3544 end do 3545 ! DEBUG 3546 ! write(std_out,*)'gmet' 3547 ! do ii=1,3 3548 ! write(std_out,*)gmet(ii,:) 3549 ! end do 3550 ! write(std_out,*)'gmet_str' 3551 ! do ii=1,2 3552 ! write(std_out,*)dtefield%gmet_str(ii,:,idir) 3553 ! end do 3554 ! ENDDEBUG 3555 do istr = 1, dtefield%nstr(idir) 3556 do ii = 1,3 3557 if (ii<idir) dtefield%coord_str(ii,istr,idir)=dtefield%fkptns(ii,dtefield%idxkstr(1,istr,idir)) 3558 if (ii>idir) dtefield%coord_str(ii-1,istr,idir)=dtefield%fkptns(ii,dtefield%idxkstr(1,istr,idir)) 3559 end do 3560 end do 3561 3562 ! the following is very similar to getshell 3563 dist_ = 0._dp 3564 do ii = 1,2 3565 dist_ = dist_ + dtefield%gmet_str(ii,ii,idir) 3566 end do 3567 max_dist = 2._dp * dist_ * 2._dp 3568 3569 dk_str(:,:,idir) = 0._dp 3570 last_dist = 0._dp 3571 ! ishell = 0 3572 ! dtefield%str_neigh(:,:,:) = 0 3573 dk_flag = 0 3574 do while (dk_flag /= 2) 3575 ! Advance shell counter 3576 ! ishell = ishell + 1 3577 3578 ! Search the smallest distance between two strings 3579 dist = max_dist 3580 do istr = 1,dtefield%nstr(idir) 3581 delta_str3(:) = dtefield%coord_str(:,1,idir) - dtefield%coord_str(:,istr,idir) 3582 do coord1 = -1,1 !two loop to search also on the border of the BZ 3583 do coord2 = -1,1 3584 dist_ = 0._dp 3585 dstr(:) = delta_str3(:) - nint(delta_str3(:)) 3586 dstr(1) = dstr(1) + real(coord1,dp) 3587 dstr(2) = dstr(2) + real(coord2,dp) 3588 do ii = 1,2 3589 do jj = 1,2 3590 dist_ = dist_ + dstr(ii)*dtefield%gmet_str(ii,jj,idir)*dstr(jj) 3591 end do 3592 end do 3593 if ((dist_ < dist).and.(dist_ - last_dist > tol8)) then 3594 dist = dist_ 3595 end if 3596 end do 3597 end do 3598 end do 3599 3600 last_dist = dist 3601 3602 ! search the connecting vectors for that distance 3603 do istr = 1,dtefield%nstr(idir) 3604 delta_str3(:) = dtefield%coord_str(:,istr,idir) - dtefield%coord_str(:,1,idir) 3605 do coord1 = -1,1 3606 do coord2 = -1,1 3607 dist_ = 0._dp 3608 dstr(:) = delta_str3(:) - nint(delta_str3(:)) 3609 dstr(1) = dstr(1) + real(coord1,dp) 3610 dstr(2) = dstr(2) + real(coord2,dp) 3611 do ii = 1,2 3612 do jj = 1,2 3613 dist_ = dist_ + dstr(ii)*dtefield%gmet_str(ii,jj,idir)*dstr(jj) 3614 end do 3615 end do 3616 if (abs(dist_ - dist) < tol8) then 3617 if (dk_flag == 0) then 3618 dk_str(:,1,idir) = dstr(:) 3619 dk_flag = 1 3620 ! DEBUG 3621 ! write(std_out,'(a,i4,2e15.4)')'1st connect', istr, dstr 3622 ! ENDDEBUG 3623 elseif (dk_str(1,1,idir)*dstr(2)-dk_str(2,1,idir)*dstr(1) > tol8) then 3624 dk_str(:,2,idir) = dstr(:) 3625 dk_flag = 2 3626 ! DEBUG 3627 ! write(std_out,'(a,i4,2e15.4)')'2nd connect', istr, dstr 3628 ! ENDDEBUG 3629 exit 3630 end if 3631 end if 3632 end do 3633 if (dk_flag == 2) exit 3634 end do 3635 if (dk_flag == 2) exit 3636 end do 3637 3638 end do ! do while 3639 3640 ! search the two neighbours for each string 3641 do istr = 1,dtefield%nstr(idir) 3642 dtefield%str_neigh(0,istr,idir) = istr 3643 dtefield%strg_neigh(0,istr,:,idir) = 0 3644 do jstr = 1,dtefield%nstr(idir) 3645 delta_str3(:) = dtefield%coord_str(:,jstr,idir) - dtefield%coord_str(:,istr,idir) 3646 do coord1 = -1,1 3647 do coord2 = -1,1 3648 dist_ = 0._dp 3649 dstr(:) = delta_str3(:) - nint(delta_str3(:)) 3650 dstr(1) = dstr(1) + real(coord1,dp) 3651 dstr(2) = dstr(2) + real(coord2,dp) 3652 do ii = 1,2 3653 if (sum(abs(dstr(:)-dk_str(:,ii,idir)))<tol8) then 3654 dtefield%str_neigh(ii,istr,idir) = jstr 3655 dtefield%strg_neigh(ii,istr,1,idir) = coord1 3656 dtefield%strg_neigh(ii,istr,2,idir) = coord2 3657 elseif (sum(abs(dstr(:)+dk_str(:,ii,idir)))<tol8) then 3658 dtefield%str_neigh(-ii,istr,idir) = jstr 3659 dtefield%strg_neigh(-ii,istr,1,idir) = coord1 3660 dtefield%strg_neigh(-ii,istr,2,idir) = coord2 3661 end if 3662 end do 3663 end do 3664 end do 3665 end do 3666 end do 3667 3668 ! DEBUG 3669 ! 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) 3670 ! write(std_out,*)'istr, neigh1, strg(1,:), neigh2, strg(2,:),neigh-1, strg(-1,:), neigh-2, strg(-2,:)' 3671 ! do istr=1,dtefield%nstr(idir) 3672 ! write(std_out,'(13i4)')istr, & 3673 ! & dtefield%str_neigh(1,istr,idir), dtefield%strg_neigh(1,istr,:,idir),& 3674 ! & dtefield%str_neigh(2,istr,idir), dtefield%strg_neigh(2,istr,:,idir),& 3675 ! & dtefield%str_neigh(-1,istr,idir), dtefield%strg_neigh(-1,istr,:,idir),& 3676 ! & dtefield%str_neigh(-2,istr,idir), dtefield%strg_neigh(-2,istr,:,idir) 3677 ! end do 3678 ! ENDDEBUG 3679 3680 3681 end if ! rfdir(idir) == 1 3682 3683 end do ! close loop over idir 3684 3685 ABI_FREE(kpt_mark) 3686 3687 call timab(1006,2,tsec) 3688 call timab(1007,1,tsec) 3689 3690 !------------------------------------------------------------------------------ 3691 !------------ Compute PAW on-site terms if necessary -------------------------- 3692 !------------------------------------------------------------------------------ 3693 3694 if (usepaw == 1 .and. dtefield%has_expibi == 1) then 3695 ABI_MALLOC(calc_expibi,(2,natom)) 3696 do idir = 1, 3 3697 dk = dtefield%dkvecs(1:3,idir) 3698 calc_expibi = zero 3699 call expibi(calc_expibi,dk,natom,xred) 3700 dtefield%expibi(1:2,1:natom,idir) = calc_expibi 3701 end do 3702 ! call expibi(dtefield%expibi,dtefield%dkvecs,natom,xred) 3703 dtefield%has_expibi = 2 3704 ABI_FREE(calc_expibi) 3705 end if 3706 3707 if (usepaw == 1 .and. dtefield%has_qijb == 1) then 3708 ABI_MALLOC(calc_qijb,(2,dtefield%lmn2max,natom)) 3709 3710 do idir = 1, 3 3711 dk = dtefield%dkvecs(1:3,idir) 3712 calc_qijb = zero 3713 call qijb_kk(calc_qijb,dk,dtefield%expibi(1:2,1:natom,idir),& 3714 & gprimd,dtefield%lmn2max,natom,ntypat,pawang,pawrad,pawtab,typat) 3715 dtefield%qijb_kk(1:2,1:dtefield%lmn2max,1:natom,idir) = calc_qijb 3716 ! call qijb_kk(dtefield%qijb_kk,dtefield%dkvecs,dtefield%expibi,& 3717 ! & gprimd,dtefield%lmn2max,natom,ntypat,pawang,pawrad,pawtab,typat) 3718 end do 3719 dtefield%has_qijb = 2 3720 ABI_FREE(calc_qijb) 3721 end if 3722 3723 if (usepaw == 1 .and. dtefield%has_rij == 1) then 3724 c1=sqrt(four_pi/three) 3725 do itypat = 1, ntypat 3726 do klmn = 1, pawtab(itypat)%lmn2_size 3727 dtefield%rij(klmn,itypat,1) = c1*pawtab(itypat)%qijl(4,klmn) ! S_{1,1} ~ x 3728 dtefield%rij(klmn,itypat,2) = c1*pawtab(itypat)%qijl(2,klmn) ! S_{1,-1} ~ y 3729 dtefield%rij(klmn,itypat,3) = c1*pawtab(itypat)%qijl(3,klmn) ! S_{1,0} ~ z 3730 end do ! end loop over klmn 3731 end do ! end loop over itypat 3732 dtefield%has_rij = 2 3733 end if ! 3734 3735 call timab(1007,2,tsec) 3736 call timab(1008,1,tsec) 3737 3738 !------------------------------------------------------------------------------ 3739 !------------ Build the array pwind that is needed to compute the ------------- 3740 !------------ overlap matrices at k +- dk ------------- 3741 !------------------------------------------------------------------------------ 3742 3743 ecut_eff = dtset%ecut*(dtset%dilatmx)**2 3744 exchn2n3d = 0 ; istwf_k = 1 ; ikg1 = 0 3745 pwind(:,:,:) = 0 3746 pwnsfac(1,:) = 1.0_dp 3747 pwnsfac(2,:) = 0.0_dp 3748 ABI_MALLOC(kg1_k,(3,mpw)) 3749 3750 ipwnsfac = 0 3751 3752 do idir = 1, 3 3753 3754 if (dtset%rfdir(idir) == 1) then 3755 3756 dk(:) = dtefield%dkvecs(:,idir) 3757 3758 do ifor = 1, 2 3759 3760 if (ifor == 2) dk(:) = -1._dp*dk(:) 3761 3762 ! Build pwind and kgindex 3763 ! NOTE: The array kgindex is important for parallel execution. 3764 ! In case nsppol = 2, it may happent that a particular processor 3765 ! treats k-points at different spin polarizations. 3766 ! In this case, it is not possible to address the elements of 3767 ! pwind correctly without making use of the kgindex array. 3768 3769 ikg = 0 ; ikpt_loc = 0 ; isppol = 1 3770 do ikpt = 1, dtefield%fnkpt 3771 3772 ikpti = dtefield%indkk_f2ibz(ikpt,1) 3773 nband_k = dtset%nband(ikpti) 3774 ikpt1f = dtefield%ikpt_dk(ikpt,ifor,idir) 3775 ikpt1i = dtefield%indkk_f2ibz(ikpt1f,1) 3776 3777 if ((proc_distrb_cycle(mpi_enreg%proc_distrb,ikpti,1,nband_k,1,me)).and.& 3778 & (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpti,1,nband_k,nsppol,me))) cycle 3779 3780 ikpt_loc = ikpt_loc + 1 3781 3782 ! Build basis sphere of plane waves for the nearest neighbour of 3783 ! the k-point (important for MPI //) 3784 3785 kg1_k(:,:) = 0 3786 kpt1(:) = dtset%kptns(:,ikpt1i) 3787 call kpgsph(ecut_eff,exchn2n3d,gmet,ikg1,ikpt,istwf_k,kg1_k,kpt1,& 3788 & 1,mpi_enreg,mpw,npw_k1) 3789 me_g0=mpi_enreg%me_g0 3790 3791 3792 ! ji: fkgindex is defined here ! 3793 dtefield%fkgindex(ikpt) = ikg 3794 3795 ! 3796 ! Deal with symmetry transformations 3797 ! 3798 3799 ! bra k-point k(b) and IBZ k-point kIBZ(b) related by 3800 ! k(b) = alpha(b) S(b)^t kIBZ(b) + G(b) 3801 ! where alpha(b), S(b) and G(b) are given by indkk_f2ibz 3802 ! 3803 ! For the ket k-point: 3804 ! k(k) = alpha(k) S(k)^t kIBZ(k) + G(k) - GBZ(k) 3805 ! where GBZ(k) takes k(k) to the BZ 3806 ! 3807 3808 isym = dtefield%indkk_f2ibz(ikpt,2) 3809 isym1 = dtefield%indkk_f2ibz(ikpt1f,2) 3810 3811 ! Construct transformed G vector that enters the matching condition: 3812 ! alpha(k) S(k)^{t,-1} ( -G(b) - GBZ(k) + G(k) ) 3813 3814 dg(:) = -dtefield%indkk_f2ibz(ikpt,3:5) & 3815 & -nint(-dtefield%fkptns(:,ikpt) - dk(:) - tol10 + & 3816 & dtefield%fkptns(:,ikpt1f)) & 3817 & +dtefield%indkk_f2ibz(ikpt1f,3:5) 3818 3819 ! old code 3820 ! iadum(:)=0 3821 ! do idum=1,3 3822 ! iadum(:)=iadum(:)+ symrec(:,idum,isym1)*dg(idum) 3823 ! end do 3824 3825 ! new code 3826 iadum(:) = MATMUL(TRANSPOSE(dtset%symrel(:,:,isym1)),dg(:)) 3827 3828 dg(:) = iadum(:) 3829 3830 if ( dtefield%indkk_f2ibz(ikpt1f,6) == 1 ) dg(:) = -dg(:) 3831 3832 ! Construct S(k)^{t,-1} S(b)^{t} 3833 3834 dum33(:,:) = MATMUL(TRANSPOSE(dtset%symrel(:,:,isym1)),symrec(:,:,isym)) 3835 3836 ! Construct alpha(k) alpha(b) 3837 3838 if (dtefield%indkk_f2ibz(ikpt,6) == dtefield%indkk_f2ibz(ikpt1f,6)) then 3839 itrs=0 3840 else 3841 itrs=1 3842 end if 3843 3844 3845 npw_k = npwarr(ikpti) 3846 ! npw_k1 = npwarr(ikpt1i) 3847 3848 ! loop over bra G vectors 3849 do ipw = 1, npw_k 3850 3851 ! NOTE: the bra G vector is taken for the sym-related IBZ k point, 3852 ! not for the FBZ k point 3853 iadum(:) = kg(:,dtefield%kgindex(ikpti) + ipw) 3854 3855 ! Store non-symmorphic operation phase factor exp[i2\pi \alpha G \cdot t] 3856 3857 if ( ipwnsfac == 0 ) then 3858 ! old code 3859 rdum=0.0_dp 3860 do idum=1,3 3861 rdum=rdum+dble(iadum(idum))*dtset%tnons(idum,isym) 3862 end do 3863 rdum=two_pi*rdum 3864 if ( dtefield%indkk_f2ibz(ikpt,6) == 1 ) rdum=-rdum 3865 pwnsfac(1,ikg+ipw) = cos(rdum) 3866 pwnsfac(2,ikg+ipw) = sin(rdum) 3867 ! 3868 ! new code 3869 ! rdum = DOT_PRODUCT(dble(iadum(:)),dtset%tnons(:,isym)) 3870 ! rdum= two_pi*rdum 3871 ! if ( dtefield%indkk_f2ibz(ikpt,6) == 1 ) rdum=-rdum 3872 ! pwnsfac(1,ikg+ipw) = cos(rdum) 3873 ! pwnsfac(2,ikg+ipw) = sin(rdum) 3874 3875 end if 3876 3877 ! to determine r.l.v. matchings, we transformed the bra vector 3878 ! Rotation 3879 iadum1(:)=0 3880 do idum1=1,3 3881 iadum1(:)=iadum1(:)+dum33(:,idum1)*iadum(idum1) 3882 end do 3883 iadum(:)=iadum1(:) 3884 ! Time reversal 3885 if (itrs==1) iadum(:)=-iadum(:) 3886 ! Translation 3887 iadum(:) = iadum(:) + dg(:) 3888 3889 do jpw = 1, npw_k1 3890 iadum1(1:3) = kg1_k(1:3,jpw) 3891 if ( (iadum(1) == iadum1(1)).and. & 3892 & (iadum(2) == iadum1(2)).and. & 3893 & (iadum(3) == iadum1(3)) ) then 3894 pwind(ikg + ipw,ifor,idir) = jpw 3895 ! write(std_out,'(a,2x,3i4,2x,i4)') 'Found !:',iadum1(:),jpw 3896 exit 3897 end if 3898 end do 3899 end do 3900 3901 ikg = ikg + npw_k 3902 3903 end do ! close loop over ikpt 3904 3905 ipwnsfac = 1 3906 3907 end do ! close loop over ifor 3908 3909 end if ! rfdir(idir) == 1 3910 3911 end do ! close loop over idir 3912 3913 3914 call timab(1008,2,tsec) 3915 call timab(1009,1,tsec) 3916 3917 !Build mpi_enreg%kptdstrb 3918 !array required to communicate the WFs between cpus in berryphase_new.f 3919 !(MPI // over k-points) 3920 if (nproc>1) then 3921 do idir = 1, 3 3922 if (dtset%rfdir(idir) == 1) then 3923 do ifor = 1, 2 3924 3925 ikpt_loc = 0 3926 do isppol = 1, nsppol 3927 3928 do ikpt = 1, dtefield%fnkpt 3929 3930 ikpti = dtefield%indkk_f2ibz(ikpt,1) 3931 nband_k = dtset%nband(ikpti) 3932 ikpt1f = dtefield%ikpt_dk(ikpt,ifor,idir) 3933 ikpt1i = dtefield%indkk_f2ibz(ikpt1f,1) 3934 3935 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpti,1,nband_k,isppol,me)) cycle 3936 3937 ikpt_loc = ikpt_loc + 1 3938 mpi_enreg%kptdstrb(me + 1,ifor+2*(idir-1),ikpt_loc) = & 3939 & ikpt1i + (isppol - 1)*nkpt 3940 3941 mpi_enreg%kptdstrb(me+1,ifor+2*(idir-1),& 3942 & ikpt_loc+dtefield%fmkmem_max*nsppol) = & 3943 & ikpt1f + (isppol - 1)*dtefield%fnkpt 3944 3945 end do ! ikpt 3946 end do ! isppol 3947 end do ! ifor 3948 end if ! dtset%rfdir(idir) == 1 3949 end do ! idir 3950 end if ! nproc>1 3951 3952 !build mpi_enreg%kpt_loc2fbz_sp 3953 ikpt_loc = 0 3954 do isppol = 1, nsppol 3955 do ikpt = 1, dtefield%fnkpt 3956 3957 ikpti = dtefield%indkk_f2ibz(ikpt,1) 3958 nband_k = dtset%nband(ikpti) 3959 3960 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpti,1,nband_k,isppol,me)) cycle 3961 3962 ikpt_loc = ikpt_loc + 1 3963 3964 mpi_enreg%kpt_loc2fbz_sp(me, ikpt_loc, 1) = ikpt 3965 mpi_enreg%kpt_loc2fbz_sp(me, ikpt_loc, 2) = isppol 3966 3967 end do 3968 end do 3969 3970 3971 !parallel case only : 3972 !build mpi_enreg%kpt_loc2ibz_sp, dtefield%cgqindex and dtefield%nneigh 3973 if ((fieldflag).and.(nproc>1)) then 3974 ikpt_loc = 0 3975 do isppol = 1, nsppol 3976 do ikpt = 1, nkpt 3977 3978 ikptf = dtefield%i2fbz(ikpt) 3979 nband_k = dtset%nband(ikpti) 3980 3981 neigh(:) = 0 ; icg = 0 ; ikg = 0 ; flag_kpt = 0; icprj = 0 3982 do idir=1, 3 3983 3984 ! skip idir values for which efield_dot(idir) = 0 3985 if (abs(dtefield%efield_dot(idir)) < tol12) cycle 3986 3987 do ifor = 1, 2 3988 3989 flag = 0 3990 3991 ikpt1f = dtefield%ikpt_dk(ikptf,ifor,idir) 3992 ikpt1i = dtefield%indkk_f2ibz(ikpt1f,1) 3993 3994 dtefield%cgqindex(3,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = ikg 3995 ikg = ikg + npwarr(ikpt1i) 3996 3997 ! check if this neighbour is also a previous neighbour 3998 do ineigh = 1, (ifor+2*(idir-1)) 3999 if (neigh(ineigh) == ikpt1i) then 4000 flag = 1 4001 dtefield%cgqindex(1,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = ineigh 4002 dtefield%cgqindex(2,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = & 4003 & dtefield%cgqindex(2,ineigh,ikpt+(isppol-1)*nkpt) 4004 exit 4005 end if 4006 end do 4007 ! create the cgqindex of the neighbour if necessary 4008 if (flag == 0) then 4009 neigh(ifor+2*(idir-1)) = ikpt1i 4010 dtefield%cgqindex(1,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = & 4011 & ifor+2*(idir-1) 4012 dtefield%cgqindex(2,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = icg 4013 if (isppol == 1) dtefield%nneigh(ikpt) = dtefield%nneigh(ikpt) + 1 4014 icg = icg + npwarr(ikpt1i)*dtefield%nspinor*nband_k 4015 end if 4016 end do !ifor 4017 end do !idir 4018 4019 if (.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me))) then 4020 ! ikpt is one of my kpt_loc 4021 ikpt_loc = ikpt_loc + 1 4022 mpi_enreg%kpt_loc2ibz_sp(me, ikpt_loc, 1) = ikpt 4023 mpi_enreg%kpt_loc2ibz_sp(me, ikpt_loc, 2) = isppol 4024 end if 4025 4026 end do !ikpt 4027 end do !isppol 4028 end if !nproc>1 4029 4030 !should be temporary 4031 !unassigned mpi_enreg%kpt_loc2fbz_sp are empty ; inform other cpu (there are better ways...) 4032 mpi_enreg%mkmem(me) = mkmem 4033 !do ii=ikpt_loc+1,dtefield%fmkmem_max 4034 !mpi_enreg%kpt_loc2fbz_sp(me, ii, 1) = -1 4035 !end do 4036 4037 4038 !(same as mpi_enreg%kptdstrb but for k-points in the iBZ), 4039 !dtefield%cgqindex and dtefield%nneigh 4040 4041 if ((fieldflag).and.(nproc>1)) then 4042 4043 ikpt_loc = 1 4044 do isppol = 1, nsppol 4045 do ikpt = 1, nkpt 4046 4047 nband_k = dtset%nband(ikpt) 4048 ikptf = dtefield%i2fbz(ikpt) 4049 4050 neigh(:) = 0 ; icg = 0 ; ikg = 0 ; flag_kpt = 0; icprj = 0 4051 do idir = 1, 3 4052 4053 ! Skip idir values for which efield_dot(idir) = 0 4054 if (abs(dtefield%efield_dot(idir)) < tol12 .and. (fieldflag)) cycle 4055 4056 do ifor = 1, 2 4057 4058 ikpt1f = dtefield%ikpt_dk(ikptf,ifor,idir) 4059 ikpt1i = dtefield%indkk_f2ibz(ikpt1f,1) 4060 4061 ! dtefield%cgqindex(3,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = ikg 4062 ikg = ikg + npwarr(ikpt1i) 4063 4064 flag = 0 4065 do ineigh = 1, (ifor+2*(idir-1)) 4066 if (neigh(ineigh) == ikpt1i) then 4067 flag = 1 4068 ! dtefield%cgqindex(1,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = ineigh 4069 ! dtefield%cgqindex(2,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = & 4070 ! & dtefield%cgqindex(2,ineigh,ikpt+(isppol-1)*nkpt) 4071 exit 4072 end if 4073 end do 4074 if (flag == 0) then 4075 ! neigh(ifor+2*(idir-1)) = ikpt1i 4076 ! dtefield%cgqindex(1,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = & 4077 ! & ifor+2*(idir-1) 4078 ! dtefield%cgqindex(2,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = icg 4079 ! if (isppol == 1) dtefield%nneigh(ikpt) = dtefield%nneigh(ikpt) + 1 4080 ! icg = icg + npwarr(ikpt1i)*dtset%nspinor*nband_k 4081 end if 4082 4083 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) cycle 4084 4085 flag_kpt = 1 4086 4087 ! MVeithen: the if condition allows to avoid that the same wavefunction 4088 ! is send several times to a particular cpu 4089 4090 end do ! ifor 4091 end do ! idir 4092 4093 if (flag_kpt == 1) ikpt_loc = ikpt_loc + 1 4094 4095 end do ! ikpt 4096 end do ! isppol 4097 4098 end if ! fieldflag 4099 4100 call xmpi_sum(mpi_enreg%kptdstrb,spaceComm,ierr) 4101 call xmpi_sum(mpi_enreg%kpt_loc2fbz_sp,spaceComm,ierr) 4102 if (fieldflag) then 4103 call xmpi_sum(mpi_enreg%kpt_loc2ibz_sp,spaceComm,ierr) 4104 call xmpi_sum(mpi_enreg%mkmem,spaceComm,ierr) 4105 end if 4106 4107 !------------------------------------------------------------------------------ 4108 !------------------------ Estimate critical field ----------------------------- 4109 !------------------------------------------------------------------------------ 4110 4111 !Compute the minimal value of the bandgap required to be below 4112 !the critical field as defined by the relation 4113 !| E_i*a_i | < E_g/n_i 4114 4115 if (fieldflag) then 4116 4117 do idir = 1, 3 4118 ! eg_dir(idir) = abs(dtefield%efield_dot(idir))*dtefield%nkstr(idir) 4119 eg_dir(idir) = abs(dtset%red_efieldbar(idir))*dtefield%nkstr(idir) 4120 end do 4121 4122 4123 eg = maxval(eg_dir) 4124 eg_ev = eg*Ha_eV 4125 4126 if (dtset%optcell ==0 .and. (dtset%berryopt == 4 .or. dtset%berryopt == 14)) then 4127 write(message,'(a,a,a,a,a,a,a,a,f7.2,a,a)')ch10,& 4128 & ' initberry: COMMENT - ',ch10,& 4129 & ' As a rough estimate,',ch10,& 4130 & ' to be below the critical field, the bandgap of your system',ch10,& 4131 & ' should be larger than ',eg_ev,' eV.',ch10 4132 call wrtout(ab_out,message,'COLL') 4133 call wrtout(std_out,message,'COLL') 4134 4135 else 4136 4137 write(message,'(a,a,a,a,a,a,a)') ch10,& 4138 & ' initberry: COMMENT - ',ch10,& 4139 & ' The estimation of critical electric field should be checked after calculation.',ch10,& 4140 & ' It is printed out just after total energy.' ,ch10 4141 4142 call wrtout(ab_out,message,'COLL') 4143 call wrtout(std_out,message,'COLL') 4144 4145 end if 4146 4147 end if 4148 4149 ABI_FREE(kg1_k) 4150 4151 call timab(1009,2,tsec) 4152 call timab(1001,2,tsec) 4153 4154 DBG_EXIT("COLL") 4155 4156 end subroutine initberry
ABINIT/m_berryphase_new [ Modules ]
NAME
m_berryphase_new
FUNCTION
COPYRIGHT
Copyright (C) 2003-2024 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 .
SOURCE
15 #if defined HAVE_CONFIG_H 16 #include "config.h" 17 #endif 18 19 #include "abi_common.h" 20 21 module m_berryphase_new 22 23 use defs_basis 24 use defs_wvltypes 25 use m_efield 26 use m_errors 27 use m_abicore 28 use m_xmpi 29 use m_hdr 30 use m_dtset 31 use m_dtfil 32 33 use defs_datatypes, only : pseudopotential_type 34 use defs_abitypes, only : MPI_type 35 use m_berrytk, only : smatrix, polcart 36 use m_cgprj, only : ctocprj 37 use m_fftcore, only : kpgsph 38 use m_geometry, only : xred2xcart, metric 39 use m_io_tools, only : open_file 40 use m_iowf, only : outwf, outresid 41 use m_kg, only : getph 42 use m_kpts, only : listkk, smpbz 43 use m_mpinfo, only : proc_distrb_cycle 44 use m_numeric_tools,only : rhophi 45 use m_pawang, only : pawang_type 46 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_get, pawcprj_mpi_allgather, & 47 pawcprj_put, pawcprj_copy, pawcprj_mpi_recv, & 48 pawcprj_mpi_send, pawcprj_free, pawcprj_getdim, pawcprj_symkn 49 use m_paw_dfpt, only : dsdr_k_paw 50 use m_paw_efield, only : pawpolev 51 use m_pawrad, only : pawrad_type 52 use m_pawrhoij, only : pawrhoij_type 53 use m_paw_sphharm, only : setsym_ylm 54 use m_pawtab, only : pawtab_type 55 use m_paw_overlap, only : expibi,qijb_kk,smatrix_k_paw 56 use m_symtk, only : symatm 57 use m_time, only : timab 58 59 implicit none 60 61 private
ABINIT/prtefield [ Functions ]
NAME
prtefield
FUNCTION
Print components of electric field, displacement field and polarization in nice format
COPYRIGHT
Copyright (C) 1998-2024 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)
SOURCE
2318 subroutine prtefield(dtset,dtefield,iunit,rprimd) 2319 2320 !Arguments ------------------------------------ 2321 integer :: iunit 2322 real(dp),intent(in) :: rprimd(3,3) 2323 type(efield_type),intent(in) :: dtefield 2324 type(dataset_type),intent(inout) :: dtset 2325 2326 !Local variables------------------------------- 2327 ! Do not modify the length of this string 2328 !scalars 2329 integer :: idir,ii 2330 character(len=1500) :: message 2331 character(len=7) :: flag_field(3) 2332 2333 real(dp) :: ucvol 2334 ! arrays 2335 real(dp) :: ptot_cart(3),gmet(3,3),gprimd(3,3),rmet(3,3),red_pbar(3),red_dbar(3),red_dfieldbar(3) 2336 real(dp) :: red_efieldbar_lc(3),red_efield_lc(3) 2337 2338 2339 ! ************************************************************************* 2340 2341 !DEBUG 2342 !write(iout, '(a)') ' prtefield : enter ' 2343 !ENDDEBUG 2344 2345 !write here 2346 2347 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 2348 2349 2350 ptot_cart(:)=zero 2351 do idir = 1,3 2352 ptot_cart(idir)=rprimd(idir,1)*dtefield%red_ptot1(1) + rprimd(idir,2)*dtefield%red_ptot1(2) + & 2353 & rprimd(idir,3)*dtefield%red_ptot1(3) 2354 end do 2355 ptot_cart(:)=ptot_cart(:)/ucvol 2356 2357 if (dtset%berryopt == 4) then 2358 2359 ! to calculate e Eq.(25) 2360 do idir=1,3 2361 dtset%red_efield(idir) =(ucvol/(4*pi))*dot_product(dtset%efield(:),gprimd(:,idir)) 2362 end do 2363 2364 ! to calculate ebar Eq.(25) 2365 do idir=1,3 2366 dtset%red_efieldbar(idir) =dot_product(dtset%efield(:),rprimd(:,idir)) 2367 end do 2368 2369 2370 ! to calculate pbar 2371 do idir=1,3 2372 red_pbar(idir) = (4*pi/ucvol)*dot_product(dtefield%red_ptot1(:),rmet(:,idir)) 2373 end do 2374 2375 !MGNAG: This message is too long and causes 2376 ! Runtime Error: wrtout_cpp.f90, line 893: Buffer overflow on output 2377 ! with NAG in test seq_tsv6_125 where we write to std_out! 2378 ! I cannot change the RECLEN of std_out! 2379 2380 write(message,'(a,a,a,3(es16.9,2x),a,a,3(es16.9,2x))') ' (a. u.)', ch10,& 2381 & ' E: ', (dtset%efield(ii), ii=1,3), ch10, & 2382 & ' P: ', (ptot_cart(ii), ii=1,3) 2383 call wrtout(iunit,message,'COLL') 2384 2385 write(message,'(a,a,3(es16.9,2x),a,a,3(es16.9,2x))') ch10,& 2386 & ' ebar: ', (dtset%red_efieldbar(ii),ii=1,3), ch10, & !!HONG need to change 2387 & ' pbar: ', (red_pbar(ii),ii=1,3) 2388 call wrtout(iunit,message,'COLL') 2389 2390 write(message,'(a,a,3(es16.9,2x),a,a,3(es16.9,2x))') ch10,& 2391 & ' e: ', (dtset%red_efield(ii),ii=1,3), ch10, & 2392 & ' p: ', (dtefield%red_ptot1(ii), ii=1,3) 2393 call wrtout(iunit,message,'COLL') 2394 2395 write(message,'(a,a,a,a,3(es16.9,2x),a,a,3(es16.9,2x),a)') ch10,& 2396 & ' (S.I.), that is V/m for E, and C/m^2 for P', ch10, & 2397 & '- 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 2398 & ' P: ', (ptot_cart(ii)*(e_Cb)/(Bohr_Ang*1d-10)**2, ii=1,3),ch10 2399 call wrtout(iunit,message,'COLL') 2400 2401 end if ! berryopt ==4 2402 2403 2404 if (dtset%berryopt == 6) then 2405 do idir=1,3 2406 dtset%red_efield(idir) =(ucvol/(4*pi))*dot_product(dtset%efield(:),gprimd(:,idir)) 2407 end do 2408 2409 ! to calculate ebar !! Need to be changed 2410 do idir=1,3 2411 dtset%red_efieldbar(idir) = dot_product(dtset%efield(:),rprimd(:,idir)) 2412 end do 2413 2414 ! to calculate red_pbar 2415 do idir=1,3 2416 red_pbar(idir) = (4*pi/ucvol)*dot_product(dtefield%red_ptot1(:),rmet(:,idir)) 2417 end do 2418 2419 ! to calculate red_dbar 2420 do idir=1,3 2421 red_dbar(idir) = dot_product(dtset%dfield(:),rprimd(:,idir)) 2422 end do 2423 2424 ! to calculate d 2425 do idir=1,3 2426 dtset%red_dfield(idir) =(ucvol/(4*pi))*dot_product(dtset%dfield(:),gprimd(:,idir)) 2427 end do 2428 2429 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,& 2430 & ' e: ', (dtset%red_efield(ii),ii=1,3), ch10, & 2431 & ' p: ', (dtefield%red_ptot1(ii), ii=1,3), ch10, & 2432 & ' d: ', (dtset%red_dfield(ii),ii = 1, 3), ch10, & 2433 & ' e + p: ', (dtset%red_efield(ii)+dtefield%red_ptot1(ii),ii=1,3) 2434 call wrtout(iunit,message,'COLL') 2435 2436 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,& 2437 & ' ebar: ', (dtset%red_efieldbar(ii),ii=1,3), ch10, & !!HONG need to change 2438 & ' pbar: ', (red_pbar(ii),ii=1,3), ch10, & 2439 & ' dbar: ', (red_dbar(ii),ii=1,3), ch10, & 2440 & ' eba+pba: ', (dtset%red_efieldbar(ii)+red_pbar(ii),ii=1,3) 2441 call wrtout(iunit,message,'COLL') 2442 2443 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, & 2444 & ' E: ', (dtset%efield(ii), ii=1,3), ch10, & 2445 & ' P: ', (ptot_cart(ii), ii=1,3), ch10, & 2446 & ' D: ', (dtset%dfield(ii),ii = 1, 3), ch10, & 2447 & 'E+4*pi*P: ', (dtset%efield(ii)+4.0d0*pi*ptot_cart(ii),ii=1,3) 2448 call wrtout(iunit,message,'COLL') 2449 2450 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,& 2451 & ' (S.I.), that is V/m for E, and C/m^2 for P and D', ch10, & 2452 & '- 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 2453 & ' P: ', (ptot_cart(ii)*(e_Cb)/(Bohr_Ang*1d-10)**2, ii=1,3), ch10,& 2454 & ' D: ', ((1.0d0/(4*pi))*dtset%dfield(ii)*(e_Cb)/(Bohr_Ang*1d-10)**2,ii = 1, 3),ch10,& 2455 & '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 2456 ! eps0*(Ha_J/(e_Cb*Bohr_Ang*1d-10))=8.854187817620*5.14220652*1d-1 2457 call wrtout(iunit,message,'COLL') 2458 2459 !MGNAG Runtime Error: wrtout_cpp.f90, line 896: Buffer overflow on output 2460 2461 end if ! berryopt ==6 2462 2463 2464 if (dtset%berryopt == 14) then 2465 2466 do idir=1,3 ! ebar local 2467 red_efieldbar_lc(idir)=dot_product(dtefield%efield2(:),rprimd(:,idir)) 2468 end do 2469 2470 2471 do idir=1,3 2472 dtset%red_efield(idir) =(ucvol/(4*pi))*dot_product(dtefield%efield2(:),gprimd(:,idir)) 2473 end do 2474 2475 ! to calculate pbar 2476 do idir=1,3 2477 red_pbar(idir) = (4*pi/ucvol)*dot_product(dtefield%red_ptot1(:),rmet(:,idir)) 2478 end do 2479 2480 write(message,'(a,a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x))') ' (a. u.)', ch10,& 2481 & ' ebar0: ', (dtset%red_efieldbar(ii),ii=1,3), ch10, & 2482 & ' ebar: ', (red_efieldbar_lc(ii),ii=1,3), ch10, & 2483 & ' pbar: ', (red_pbar(ii),ii=1,3) 2484 call wrtout(iunit,message,'COLL') 2485 2486 write(message,'(a,a,3(es16.9,2x),a,a,3(es16.9,2x))') ch10, & 2487 & ' e: ', (dtset%red_efield(ii),ii=1,3), ch10, & 2488 & ' p: ', (dtefield%red_ptot1(ii), ii=1,3) 2489 call wrtout(iunit,message,'COLL') 2490 2491 write(message,'(a,a,3(es16.9,2x),a,a,3(es16.9,2x))') ch10, & 2492 & ' E: ', (dtefield%efield2(ii), ii=1,3), ch10, & 2493 & ' P: ', (ptot_cart(ii), ii=1,3) 2494 call wrtout(iunit,message,'COLL') 2495 2496 write(message,'(a,a,a,a,3(es16.9,2x),a,a,3(es16.9,2x),a)') ch10, & 2497 & ' (S.I.), that is V/m for E, and C/m^2 for P', ch10, & 2498 & '- 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 2499 & ' P: ', (ptot_cart(ii)*(e_Cb)/(Bohr_Ang*1d-10)**2, ii=1,3),ch10 2500 call wrtout(iunit,message,'COLL') 2501 2502 2503 end if ! berryopt ==14 2504 2505 2506 if (dtset%berryopt == 16) then 2507 2508 ! to calculate e Eq.(25) 2509 do idir=1,3 2510 dtset%red_efield(idir) =(ucvol/(4*pi))*dot_product(dtset%efield(:),gprimd(:,idir)) 2511 end do 2512 2513 ! to calculate ebar 2514 do idir=1,3 2515 dtset%red_efieldbar(idir) = dot_product(dtset%efield(:),rprimd(:,idir)) 2516 end do 2517 2518 ! to calculate pbar 2519 do idir=1,3 2520 red_pbar(idir) = (4*pi/ucvol)*dot_product(dtefield%red_ptot1(:),rmet(:,idir)) 2521 end do 2522 2523 ! to calculate dbar 2524 do idir=1,3 2525 red_dfieldbar(idir) = (4*pi/ucvol)*dot_product(dtset%red_dfield(:),rmet(:,idir)) 2526 end do 2527 2528 ! to calculate D 2529 do idir=1,3 2530 dtset%dfield(idir) =(4*pi/ucvol)*dot_product(dtset%red_dfield(:),rprimd(:,idir)) 2531 end do 2532 2533 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,& 2534 & ' e: ', (dtset%red_efield(ii),ii=1,3), ch10, & 2535 & ' p: ', (dtefield%red_ptot1(ii), ii=1,3), ch10, & 2536 & ' d: ', (dtset%red_dfield(ii),ii = 1, 3), ch10, & 2537 & ' e + p: ', (dtset%red_efield(ii)+dtefield%red_ptot1(ii),ii=1,3) 2538 call wrtout(iunit,message,'COLL') 2539 2540 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, & 2541 & ' ebar: ', (dtset%red_efieldbar(ii),ii=1,3), ch10, & 2542 & ' pbar: ', (red_pbar(ii),ii=1,3), ch10, & 2543 & ' dbar: ', (red_dfieldbar(ii),ii=1,3), ch10, & 2544 & ' eba+pba: ', (dtset%red_efieldbar(ii)+red_pbar(ii),ii=1,3) 2545 call wrtout(iunit,message,'COLL') 2546 2547 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, & 2548 & ' E: ', (dtset%efield(ii), ii=1,3), ch10, & 2549 & ' P: ', (ptot_cart(ii), ii=1,3), ch10, & 2550 & ' D: ', (dtset%dfield(ii),ii = 1, 3), ch10, & 2551 & 'E+4*pi*P: ', (dtset%efield(ii)+4.0d0*pi*ptot_cart(ii),ii=1,3) 2552 call wrtout(iunit,message,'COLL') 2553 2554 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, & 2555 & ' (S.I.), that is V/m for E, and C/m^2 for P and D', ch10, & 2556 & '- 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 2557 & ' P: ', (ptot_cart(ii)*(e_Cb)/(Bohr_Ang*1d-10)**2, ii=1,3), ch10, & 2558 & ' D: ', ((1.0d0/(4*pi))*dtset%dfield(ii)*(e_Cb)/(Bohr_Ang*1d-10)**2,ii = 1, 3), ch10, & 2559 & '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 2560 call wrtout(iunit,message,'COLL') 2561 2562 end if ! berryopt ==16 2563 2564 if (dtset%berryopt == 17) then 2565 2566 do idir=1,3 ! ebar local 2567 red_efieldbar_lc(idir)=dot_product(dtefield%efield2(:),rprimd(:,idir)) 2568 end do 2569 2570 do idir=1,3 2571 dtset%red_efield(idir) =(ucvol/(4*pi))*dot_product(dtefield%efield2(:),gprimd(:,idir)) 2572 end do 2573 2574 ! to calculate pbar 2575 do idir=1,3 2576 red_pbar(idir) = (4*pi/ucvol)*dot_product(dtefield%red_ptot1(:),rmet(:,idir)) 2577 end do 2578 2579 2580 ! do idir=1,3 2581 ! if (dtset%rfdir(idir)==1) then ! fixed ebar 2582 ! red_efieldbar_lc(idir)=dot_product(dtefield%efield2(:),rprimd(:,idir)) ! local efieldbar 2583 ! dtset%red_efield(idir) =(ucvol/(4*pi))*dot_product(dtefield%efield2(:),gprimd(:,idir)) 2584 ! red_pbar(idir) = (4*pi/ucvol)*dot_product(dtefield%red_ptot1(:),rmet(:,idir)) 2585 ! dtset%dfield(idir)=dtefield%efield2(idir)+4*pi*ptot_cart(idir) 2586 ! dtset%red_dfield(idir)=dtset%red_efield+dtefield%red_ptot1(idir) 2587 ! dtset%red_dfieldbar(idir)=red_efieldbar_lc(idir)+red_pbar(idir) 2588 ! E_lc(idir)=dtefield%efield2(idir) 2589 ! e_lc(idir)=red_efieldbar_lc(idir) 2590 ! ebar_lc(idir)=dtset%red_efieldbar(idir) 2591 ! else if (dtset%rfdir(idir)==2) then ! fixed d 2592 ! dtset%red_efield(idir) =(ucvol/(4*pi))*dot_product(dtefield%efield2(:),gprimd(:,idir)) 2593 ! dtset%red_efieldbar(idir) = dot_product(dtefield%efield2(:),rprimd(:,idir)) 2594 ! red_pbar(idir) = (4*pi/ucvol)*dot_product(dtefield%red_ptot1(:),rmet(:,idir)) 2595 ! red_dfieldbar(idir) = (4*pi/ucvol)*dot_product(dtset%red_dfield(:),rmet(:,idir)) 2596 ! dtset%dfield(idir) =(4*pi/ucvol)*dot_product(dtset%red_dfield(:),rprimd(:,idir)) 2597 ! E_lc(idir)=dtefield%efield2(idir) 2598 ! e_lc(idir)=dtset%red_efield(idir) 2599 ! ebar_lc(idir)=dtset%red_efieldbar(idir) 2600 ! end if 2601 ! enddo 2602 2603 2604 2605 do idir=1,3 2606 red_efield_lc(idir)= (ucvol/(4*pi))*dot_product(dtefield%efield2(:),gprimd(:,idir)) 2607 red_efieldbar_lc(idir)=dot_product(dtefield%efield2(:),rprimd(:,idir)) ! local efieldbar 2608 red_pbar(idir) = (4*pi/ucvol)*dot_product(dtefield%red_ptot1(:),rmet(:,idir)) 2609 red_dfieldbar(idir) = (4*pi/ucvol)*dot_product(dtset%red_dfield(:),rmet(:,idir)) 2610 dtset%dfield(idir) =(4*pi/ucvol)*dot_product(dtset%red_dfield(:),rprimd(:,idir)) 2611 2612 end do 2613 2614 do idir=1,3 2615 if(dtset%jfielddir(idir)==1) then 2616 flag_field(idir)="E-field" 2617 else 2618 flag_field(idir)="D-field" 2619 end if 2620 end do 2621 2622 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))') & 2623 & ' (a. u.)', ch10,& 2624 & ' ', (flag_field(ii),ii=1,3),ch10, & 2625 & ' ebar0: ', (dtset%red_efieldbar(ii),ii=1,3), ch10, & 2626 & ' ebar: ', (red_efieldbar_lc(ii),ii=1,3), ch10, & 2627 & ' d: ', (dtset%red_dfield(ii),ii = 1, 3), ch10, & 2628 & ' e + p: ', (red_efield_lc(ii)+dtefield%red_ptot1(ii),ii=1,3) 2629 call wrtout(iunit,message,'COLL') 2630 2631 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, & 2632 & ' e: ', (red_efield_lc(ii),ii=1,3), ch10, & 2633 & ' p: ', (dtefield%red_ptot1(ii), ii=1,3), ch10, & 2634 & ' d: ', (dtset%red_dfield(ii),ii = 1, 3), ch10, & 2635 & ' e + p: ', (red_efield_lc(ii)+dtefield%red_ptot1(ii),ii=1,3) 2636 call wrtout(iunit,message,'COLL') 2637 2638 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, & 2639 & ' ebar: ', (red_efieldbar_lc(ii),ii=1,3), ch10, & 2640 & ' pbar: ', (red_pbar(ii),ii=1,3), ch10, & 2641 & ' dbar: ', (red_dfieldbar(ii),ii=1,3), ch10, & 2642 & ' eba+pba: ', (red_efieldbar_lc(ii)+red_pbar(ii),ii=1,3) 2643 call wrtout(iunit,message,'COLL') 2644 2645 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, & 2646 & ' E: ', (dtefield%efield2(ii), ii=1,3), ch10, & 2647 & ' P: ', (ptot_cart(ii), ii=1,3), ch10, & 2648 & ' D: ', (dtset%dfield(ii),ii = 1, 3), ch10, & 2649 & 'E+4*pi*P: ', (dtset%efield(ii)+4.0d0*pi*ptot_cart(ii),ii=1,3) 2650 call wrtout(iunit,message,'COLL') 2651 2652 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, & 2653 & ' (S.I.), that is V/m for E, and C/m^2 for P and D', ch10, & 2654 & ' 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 2655 & ' P: ', (ptot_cart(ii)*(e_Cb)/(Bohr_Ang*1d-10)**2, ii=1,3), ch10, & 2656 & ' D: ', ((1.0d0/(4*pi))*dtset%dfield(ii)*(e_Cb)/(Bohr_Ang*1d-10)**2,ii = 1, 3), ch10, & 2657 & '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 2658 call wrtout(iunit,message,'COLL') 2659 2660 end if ! berryopt ==17 2661 2662 end subroutine prtefield
ABINIT/update_e_field_vars [ Functions ]
NAME
update_e_field_vars
FUNCTION
This routine updates E field variables
COPYRIGHT
Copyright (C) 2003-2024 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
SOURCE
1807 subroutine update_e_field_vars(atindx,atindx1,cg,dimcprj,dtefield,dtfil,dtset,& 1808 & efield_old_cart,gmet,gprimd,hdr,idir,kg,mcg,& 1809 & mkmem,mpi_enreg,mpw,my_natom,natom,nattyp,ngfft,nkpt,npwarr,ntypat,& 1810 & pawrhoij,pawtab,pel_cg,pelev,pion,psps,ptot,ptot_cart,pwind,& 1811 & pwind_alloc,pwnsfac,red_efield2,red_efield2_old,red_ptot,rmet,rprimd,& 1812 & scfcv_level,scfcv_quit,scfcv_step,ucvol,unit_out,& 1813 & usepaw,xred,ylm,ylmgr) 1814 1815 !Arguments ------------------------------------ 1816 integer, intent(in) :: idir,mcg,mkmem,mpw,my_natom,natom,nkpt,ntypat 1817 integer, intent(in) :: pwind_alloc,scfcv_level,scfcv_quit,scfcv_step,unit_out,usepaw 1818 real(dp), intent(in) :: ucvol 1819 type(datafiles_type), intent(in) :: dtfil 1820 type(pseudopotential_type),intent(in) :: psps 1821 type(dataset_type), intent(inout) :: dtset 1822 type(efield_type), intent(inout) :: dtefield 1823 type(hdr_type), intent(inout) :: hdr 1824 type(MPI_type), intent(inout) :: mpi_enreg 1825 !arrays 1826 integer, intent(in) :: atindx(natom),atindx1(natom),dimcprj(usepaw*natom) 1827 integer, intent(in) :: kg(3,mpw*mkmem),nattyp(ntypat) 1828 integer, intent(in) :: ngfft(18),npwarr(nkpt),pwind(pwind_alloc,2,3) 1829 real(dp), intent(in) :: cg(2,mcg),gmet(3,3),gprimd(3,3) 1830 real(dp), intent(in) :: pwnsfac(2,pwind_alloc) 1831 real(dp), intent(in) :: rmet(3,3),rprimd(3,3) 1832 real(dp), intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm) 1833 real(dp), intent(in) :: ylmgr(dtset%mpw*dtset%mkmem,3,psps%mpsang*psps%mpsang*psps%useylm) 1834 real(dp), intent(inout) :: ptot_cart(3),xred(3,natom),efield_old_cart(3) !vz_i 1835 real(dp), intent(out) :: pel_cg(3),pelev(3),pion(3) !vz_i 1836 real(dp), intent(inout) :: red_efield2(3),red_efield2_old(3) !vz_i 1837 real(dp), intent(out) :: ptot(3),red_ptot(3) !vz_i 1838 type(pawrhoij_type), intent(in) :: pawrhoij(my_natom*usepaw) 1839 type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw) 1840 1841 !Local variables ------------------------- 1842 !scalars 1843 character(len=500) :: message 1844 integer :: ctocprj_choice,iatom,ii,iorder_cprj,mcg13,mcprj,my_nspinor,ncpgr 1845 integer :: optberry,usecprj 1846 logical :: calc_epaw3_force, calc_epaw3_stress, efield, save_cg13 1847 !arrays 1848 real(dp) :: efield_test_cart(3),red_efield1(3) 1849 real(dp),allocatable :: cg13(:,:,:),ph1d(:,:) 1850 type(pawcprj_type),allocatable :: cprj(:,:) 1851 1852 ! ************************************************************************* 1853 1854 efield = .false. 1855 1856 if ( dtset%berryopt == 4 .or. & 1857 & dtset%berryopt == 6 .or. & 1858 & dtset%berryopt == 7 .or. & 1859 & dtset%berryopt ==14 .or. & 1860 & dtset%berryopt ==16 .or. & 1861 & dtset%berryopt ==17 ) efield = .true. 1862 calc_epaw3_force = ( efield .and. dtset%optforces /= 0 .and. usepaw == 1 ) 1863 calc_epaw3_stress = ( efield .and. dtset%optstress /= 0 .and. usepaw == 1 ) 1864 1865 usecprj=1; if (psps%usepaw==0) usecprj = 0 1866 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 1867 mcprj=my_nspinor*dtset%mband*dtset%mkmem*dtset%nsppol/mpi_enreg%nproc_band 1868 !Just to be sure 1869 if (efield.and.mpi_enreg%nproc_band/=1) then 1870 message='Berry phase calculations are not designed for band parallelization!' 1871 ABI_BUG(message) 1872 end if 1873 1874 ncpgr = 0 1875 ctocprj_choice = 1 ! no derivs 1876 if ( efield .and. psps%usepaw == 1) then 1877 ABI_MALLOC(cprj,(dtset%natom,mcprj)) 1878 ! finite electric field may need gradients for forces, stress 1879 if (calc_epaw3_force .and. .not. calc_epaw3_stress) then 1880 ncpgr = 3; ctocprj_choice = 2 ! derivs w.r.t. position 1881 else if (.not. calc_epaw3_force .and. calc_epaw3_stress) then 1882 ncpgr = 6; ctocprj_choice = 3 ! derivs w.r.t strain 1883 else if (calc_epaw3_force .and. calc_epaw3_stress) then 1884 ncpgr = 9; ctocprj_choice = 23 ! derivs w.r.t. position and strain 1885 end if 1886 call pawcprj_alloc(cprj,ncpgr,dimcprj) 1887 iatom=0 ; iorder_cprj=1 ! retain ordering of input list 1888 ! all arguments to ctocprj are defined already except ph1d, do that here 1889 ABI_MALLOC(ph1d,(2,3*(2*dtset%mgfft+1)*dtset%natom)) 1890 call getph(atindx,dtset%natom,ngfft(1),ngfft(2),ngfft(3),ph1d,xred) 1891 call ctocprj(atindx,cg,ctocprj_choice,cprj,gmet,gprimd,iatom,idir,iorder_cprj,& 1892 & dtset%istwfk,kg,dtset%kptns,mcg,mcprj,dtset%mgfft,dtset%mkmem,& 1893 & mpi_enreg,psps%mpsang,dtset%mpw,dtset%natom,nattyp,dtset%nband,& 1894 & dtset%natom,ngfft,dtset%nkpt,dtset%nloalg,npwarr,dtset%nspinor,& 1895 & dtset%nsppol,dtset%nsppol,dtset%ntypat,dtset%paral_kgb,ph1d,psps,rmet,& 1896 & dtset%typat,ucvol,dtfil%unpaw,xred,ylm,ylmgr) 1897 ABI_FREE(ph1d) 1898 else 1899 ABI_MALLOC(cprj,(0,0)) 1900 end if ! end update of cprj 1901 1902 if ( efield ) then ! compute polarization and if necessary store cprj in efield 1903 optberry=1 1904 pel_cg(:) = zero;pelev=zero 1905 1906 mcg13 = 0 1907 save_cg13 = .FALSE. 1908 ABI_MALLOC(cg13,(2,mcg13,3)) 1909 call berryphase_new(atindx1,cg,cg13,cprj,dtefield,dtfil,dtset,psps,gprimd,hdr,psps%indlmn,kg,& 1910 & psps%lmnmax,dtset%mband,mcg,mcg13,mcprj,dtset%mkmem,mpi_enreg,dtset%mpw,my_natom,& 1911 & dtset%natom,npwarr,dtset%nsppol,psps%ntypat,dtset%nkpt,optberry,pawrhoij,pawtab,& 1912 & pel_cg,pelev,pion,ptot,red_ptot,pwind,& 1913 & pwind_alloc,pwnsfac,rprimd,save_cg13,dtset%typat,ucvol,& 1914 & unit_out,usecprj,psps%usepaw,xred,psps%ziontypat) 1915 ABI_FREE(cg13) 1916 1917 dtefield%red_ptot1(:)=red_ptot(:) 1918 1919 end if ! end compute polarization and store cprj for efield 1920 1921 if (efield .and. (scfcv_level == 0) ) then ! do this before scfcv loop 1922 1923 efield_old_cart(:)=dtset%efield(:) !!HONG 1924 1925 ! save this value in order to print the final value of real electric field, comparing with the desired red_fieldbar 1926 dtefield%efield2(:)=dtset%efield(:) 1927 1928 if ( dtset%berryopt ==16 .or. dtset%berryopt ==17) then !!HONG 1929 do ii=1,3 1930 red_efield2(ii)=zero 1931 red_efield2_old(ii) =(ucvol/(4*pi))*dot_product(dtset%efield(:),gprimd(:,ii)) 1932 end do 1933 end if 1934 1935 if (dtset%berryopt == 14 .and. scfcv_quit /=1) then 1936 ! ! Convert polarization to cartesian coords 1937 1938 ptot_cart(:)=zero 1939 do ii = 1,3 1940 ptot_cart(ii)=rprimd(ii,1)*red_ptot(1) + rprimd(ii,2)*red_ptot(2) + & 1941 & rprimd(ii,3)*red_ptot(3) 1942 end do 1943 ptot_cart(:)=ptot_cart(:)/ucvol 1944 1945 do ii=1,3 1946 dtefield%efield_dot(ii) = dot_product(dtset%efield(:),rprimd(:,ii)) 1947 end do 1948 1949 ! !write the field parameters: D, E, P, d, e, p, dbar, ebar, pbar 1950 write(message,'(a,a)') ch10, 'scfcv: Constant reduced ebar-field:' 1951 1952 call wrtout(std_out,message,'COLL') 1953 call prtefield(dtset,dtefield,std_out,rprimd) 1954 1955 if(dtset%prtvol>=10)then 1956 call wrtout(ab_out,message,'COLL') 1957 call prtefield(dtset,dtefield,ab_out,rprimd) 1958 end if 1959 1960 ! updating E field 1961 do ii =1,3 ! desired E field 1962 efield_test_cart(ii)=gprimd(ii,1)*dtset%red_efieldbar(1) + & 1963 & gprimd(ii,2)*dtset%red_efieldbar(2)+gprimd(ii,3)*dtset%red_efieldbar(3) 1964 end do 1965 1966 ! if not convergence well, need to add some code here to make sure efield_test_cart(:) not change much 1967 dtset%efield(:) = efield_test_cart(:) 1968 1969 end if ! berryopt ==14 1970 1971 end if ! end efield .and. scfcv_level 0 tasks 1972 1973 !!! 1974 !!! Various printing and update steps for the different efield options 1975 !!! 1976 1977 if (efield .and. (scfcv_level == 1) ) then ! do this each scf step 1978 1979 if (dtset%prtvol >= 10)then 1980 write(message,'(6(a),3(e16.9,2x),a,a,3(e16.9,2x))')ch10,& 1981 & ' scfcv: New value of the polarization:',ch10,& 1982 & ' (reduced coordinates, a. u.)',ch10,& 1983 & ' Electronic berry phase: ', (pel_cg(ii), ii = 1, 3) 1984 call wrtout(ab_out,message,'COLL') 1985 call wrtout(std_out,message,'COLL') 1986 if(psps%usepaw==1) then 1987 write(message,'(a,3(e16.9,2x))')& 1988 & ' ...includes PAW on-site term: ', (pelev(ii), ii = 1, 3) 1989 call wrtout(ab_out,message,'COLL') 1990 call wrtout(std_out,message,'COLL') 1991 end if 1992 write(message,'(a,3(e16.9,2x),a,a,3(e16.9,2x))')& 1993 & ' Ionic: ', (pion(ii), ii = 1, 3), ch10, & 1994 & ' Total: ', (red_ptot(ii), ii = 1, 3) !!REC 1995 call wrtout(ab_out,message,'COLL') 1996 call wrtout(std_out,message,'COLL') 1997 end if ! end prtvol >= 10 output 1998 1999 ptot_cart(:)=zero 2000 do ii = 1,3 2001 ptot_cart(ii)=rprimd(ii,1)*red_ptot(1) + rprimd(ii,2)*red_ptot(2) + & 2002 & rprimd(ii,3)*red_ptot(3) 2003 end do 2004 ptot_cart(:)=ptot_cart(:)/ucvol 2005 2006 ! !=================================================================================================== 2007 ! ! OUTPUT for fixed E 2008 ! !=================================================================================================== 2009 2010 if (dtset%berryopt == 4) then 2011 2012 ! !write the field parameters: D, E, P, d, e, p, dbar, ebar, pbar 2013 write(message,'(a,a)') ch10, 'scfcv: Constant unreduced E-field:' 2014 call wrtout(std_out,message,'COLL') 2015 call prtefield(dtset,dtefield,std_out,rprimd) 2016 if(dtset%prtvol>=10)then 2017 call wrtout(ab_out,message,'COLL') 2018 call prtefield(dtset,dtefield,ab_out,rprimd) 2019 end if 2020 end if ! end berryopt 4 output 2021 2022 ! ===================================================================================== 2023 ! ! fixed D calculation 2024 ! !==================================================================================== 2025 if (dtset%berryopt == 6) then 2026 if (scfcv_step > 1) then 2027 2028 ! ! update efield taking damping into account dfield is in cartesian in dtset structure (contains input value) 2029 ! ! same goes for efield - update the dtset%efield value 2030 efield_test_cart(:)=dtset%ddamp*(dtset%dfield(:)-4.0d0*pi*ptot_cart(:))+& 2031 & (1.0d0-dtset%ddamp)*efield_old_cart(:) 2032 2033 ! ! test whether change in efield in any direction exceed maxestep, if so, set the 2034 ! ! change to maxestep instead ! need optimized ! 2035 do ii = 1,3 2036 2037 if (dabs(efield_test_cart(ii)-efield_old_cart(ii)) > dabs(dtset%maxestep)) then 2038 2039 write(std_out,'(a,a,i5)') "JH - "," E-field component:",ii 2040 write(std_out,'(a,es13.5,a,es13.5,a,es13.5,a,es13.5)') " E(n)=",efield_test_cart(ii), & 2041 & ", E(n-1)=",efield_old_cart(ii), ", E(n)-E(n-1)=", efield_test_cart(ii)-efield_old_cart(ii), & 2042 & ", maxestep=",dtset%maxestep 2043 2044 2045 if (efield_test_cart(ii) > efield_old_cart(ii)) then 2046 efield_test_cart(ii) = efield_old_cart(ii) + dabs(dtset%maxestep) 2047 else 2048 efield_test_cart(ii) = efield_old_cart(ii) - dabs(dtset%maxestep) 2049 end if 2050 end if 2051 end do 2052 2053 dtset%efield(:) = efield_test_cart(:) 2054 2055 ! !write the field parameters: D, E, P, d, e, p, dbar, ebar, pbar 2056 write(message,'(a,a)') ch10, 'scfcv: Constant unreduced D-field - updating E-field:' 2057 call wrtout(std_out,message,'COLL') 2058 call prtefield(dtset,dtefield,std_out,rprimd) 2059 if(dtset%prtvol>=10)then 2060 call wrtout(ab_out,message,'COLL') 2061 call prtefield(dtset,dtefield,ab_out,rprimd) 2062 end if 2063 2064 ! ! need to update dtset%efield_dot(:) with new value 2065 dtefield%efield_dot(1) = dot_product(dtset%efield(:),rprimd(:,1)) 2066 dtefield%efield_dot(2) = dot_product(dtset%efield(:),rprimd(:,2)) 2067 dtefield%efield_dot(3) = dot_product(dtset%efield(:),rprimd(:,3)) 2068 2069 else 2070 2071 write(message,'(a,a)') ch10, 'scfcv: Constant unreduced D-field - Pre E-field:' 2072 call wrtout(std_out,message,'COLL') 2073 call prtefield(dtset,dtefield,std_out,rprimd) 2074 if(dtset%prtvol>=10)then 2075 call wrtout(ab_out,message,'COLL') 2076 call prtefield(dtset,dtefield,ab_out,rprimd) 2077 end if 2078 2079 end if ! scfcv_step >1 2080 2081 efield_old_cart(:)=dtset%efield(:) 2082 end if ! berryopt ==6 2083 ! !=================================================================================================== 2084 ! ! fixed reduced d calculation 2085 ! !=================================================================================================== 2086 if (dtset%berryopt == 16) then 2087 2088 if (scfcv_step > 1) then 2089 ! ! update efield taking damping into account reduced red_dfield 2090 ! red_efield2 is reduced electric field, defined by Eq.(25) of Nat. Phys. suppl. (2009) [[cite:Stengel2009]] 2091 2092 red_efield2(:)=dtset%ddamp*(dtset%red_dfield(:)-red_ptot(:))+ (1.0d0-dtset%ddamp)*red_efield2_old(:) 2093 2094 ! to calculate unreduced E 2095 efield_test_cart(:)=(4*pi/ucvol)*(rprimd(:,1)*red_efield2(1)+rprimd(:,2)*red_efield2(2)+rprimd(:,3)*red_efield2(3)) 2096 2097 ! ! test whether change in efield in any direction exceed maxestep, if so, set the 2098 ! ! change to maxestep instead ! need optimized ! 2099 do ii = 1,3 2100 2101 if (dabs(efield_test_cart(ii)-efield_old_cart(ii)) > dabs(dtset%maxestep)) then 2102 2103 write(std_out,'(a,a,i5)') "JH - "," E-field component:",ii 2104 write(std_out,'(a,es13.5,a,es13.5,a,es13.5,a,es13.5)') " E(n)=",efield_test_cart(ii), & 2105 & ", E(n-1)=",efield_old_cart(ii), ", E(n)-E(n-1)=", efield_test_cart(ii)-efield_old_cart(ii), & 2106 & ", maxestep=",dtset%maxestep 2107 2108 if (efield_test_cart(ii) > efield_old_cart(ii)) then 2109 efield_test_cart(ii) = efield_old_cart(ii) + dabs(dtset%maxestep) 2110 else 2111 efield_test_cart(ii) = efield_old_cart(ii) - dabs(dtset%maxestep) 2112 end if 2113 end if 2114 end do 2115 2116 dtset%efield(:) = efield_test_cart(:) 2117 2118 ! !write the field parameters: D, E, P, d, e, p, dbar, ebar, pbar 2119 write(message,'(a,a)') ch10, 'scfcv: Constant reduced d-field - updating E-field:' 2120 call wrtout(std_out,message,'COLL') 2121 call prtefield(dtset,dtefield,std_out,rprimd) 2122 if(dtset%prtvol>=10)then 2123 call wrtout(ab_out,message,'COLL') 2124 call prtefield(dtset,dtefield,ab_out,rprimd) 2125 end if 2126 2127 ! ! need to update dtset%efield_dot(:) with new value 2128 ! ! This needs to be deleted when efield_dot is deleted 2129 dtefield%efield_dot(1) = dot_product(dtset%efield(:),rprimd(:,1)) 2130 dtefield%efield_dot(2) = dot_product(dtset%efield(:),rprimd(:,2)) 2131 dtefield%efield_dot(3) = dot_product(dtset%efield(:),rprimd(:,3)) 2132 2133 else 2134 2135 write(message,'(a,a)') ch10, 'scfcv: Constant reduced d-field - Pre E-field:' 2136 call wrtout(std_out,message,'COLL') 2137 call prtefield(dtset,dtefield,std_out,rprimd) 2138 if(dtset%prtvol>=10)then 2139 call wrtout(ab_out,message,'COLL') 2140 call prtefield(dtset,dtefield,ab_out,rprimd) 2141 end if 2142 2143 end if ! scfcv_step > 1 2144 2145 efield_old_cart(:)=dtset%efield(:) 2146 red_efield2_old(:)=red_efield2(:) 2147 end if ! berryopt ==16 2148 2149 2150 ! !=================================================================================================== 2151 ! ! fixed reduced d and ebar calculation (mixed BC) 2152 ! !=================================================================================================== 2153 if (dtset%berryopt == 17) then 2154 2155 if (scfcv_step > 1) then 2156 ! ! update efield taking damping into account reduced red_dfield 2157 ! red_efield1 and red_efield2 is reduced electric field, defined by Eq.(25) of Nat. Phys. suppl. (2009) [[cite:Stengel1999]] 2158 ! red_efield1 for fixed ebar, red_efield2 for fixed d calculation 2159 2160 ! save this value in order to print the final value of real electric field, comparing with the desired red_fieldbar 2161 dtefield%efield2(:)=dtset%efield(:) 2162 2163 ! write(*,'(a,3i4)') "jfielddir=", (dtset%jfielddir(ii),ii=1,3) 2164 2165 do ii=1,3 2166 if (dtset%jfielddir(ii) ==2 ) then ! direction under fixed d 2167 dtset%red_efieldbar(ii) = dot_product(dtset%efield(:),rprimd(:,ii)) ! update ebar which is not fixed 2168 dtefield%efield_dot(ii) = dot_product(dtset%efield(:),rprimd(:,ii)) 2169 red_efield2(ii)=dtset%ddamp*(dtset%red_dfield(ii) - red_ptot(ii)) + & 2170 & (1.0d0-dtset%ddamp)*red_efield2_old(ii) ! d(ii) is fixed, update e(ii) may need ddamping here 2171 2172 ! write(message,'(a,a,i5,a,i5)') ch10, 'direction ', ii,' for fixed d, value is (2) ', dtset%jfielddir(ii) 2173 ! call wrtout(ab_out,message,'COLL') 2174 ! call wrtout(std_out,message,'COLL') 2175 2176 else if (dtset%jfielddir(ii) ==1 ) then ! direction under fixed ebar 2177 red_efield2(ii)= (ucvol/(4*pi))*dot_product(dtset%efield(:),gprimd(:,ii)) ! update e which is not fixed 2178 dtset%red_dfield(ii)=red_ptot(ii) + (ucvol/(4*pi))*dot_product(dtset%efield(:),gprimd(:,ii)) ! update d 2179 2180 ! write(message,'(a,a,i5,a,i5)') ch10, 'direction ', ii,' for fixed ebar, value is (1) ', dtset%jfielddir(ii) 2181 ! call wrtout(ab_out,message,'COLL') 2182 ! call wrtout(std_out,message,'COLL') 2183 2184 end if 2185 end do 2186 2187 do ii=1,3 2188 red_efield1(ii) =(ucvol/(4*pi))*dot_product(dtset%red_efieldbar(:),gmet(:,ii)) 2189 end do 2190 2191 2192 dtset%red_efield(:)=(red_efield1(:) + red_efield2(:))/2.0d0 ! average reduced efield, 2193 ! one is from fixed ebar part, 2194 ! the other is from fixed d part. 2195 ! This may need to be optimized !! 2196 2197 write(message,'(a,a,a,a,3(es16.9,2x),a)') ch10, 'Reduced efield from fixed ebar:', ch10, & 2198 & ' e: ', (red_efield1(ii),ii=1,3), ch10 2199 2200 ! call wrtout(ab_out,message,'COLL') 2201 call wrtout(std_out,message,'COLL') 2202 2203 write(message,'(a,a,a,a,3(es16.9,2x),a)') ch10, 'Reduced efield from fixed d:', ch10, & 2204 & ' e: ', (red_efield2(ii),ii=1,3), ch10 2205 2206 ! call wrtout(ab_out,message,'COLL') 2207 call wrtout(std_out,message,'COLL') 2208 2209 write(message,'(a,a,a,a,3(es16.9,2x),a)') ch10, 'Average reduced efield:', ch10, & 2210 & ' e: ', (dtset%red_efield(ii),ii=1,3), ch10 2211 2212 ! call wrtout(ab_out,message,'COLL') 2213 call wrtout(std_out,message,'COLL') 2214 2215 ! to calculate unreduced E 2216 do ii=1,3 2217 efield_test_cart(ii) = (4*pi/ucvol)* dot_product(dtset%red_efield(:),rprimd(:,ii)) 2218 end do 2219 2220 ! ! test whether change in efield in any direction exceed maxestep, if so, set the 2221 ! ! change to maxestep instead ! need optimized ! 2222 do ii = 1,3 2223 if (dabs(efield_test_cart(ii)-efield_old_cart(ii)) > dabs(dtset%maxestep)) then 2224 2225 write(std_out,'(a,a,i5)') "JH - "," E-field component:",ii 2226 write(std_out,'(a,es13.5,a,es13.5,a,es13.5,a,es13.5)') " E(n)=",efield_test_cart(ii), & 2227 & ", E(n-1)=",efield_old_cart(ii), ", E(n)-E(n-1)=", efield_test_cart(ii)-efield_old_cart(ii), & 2228 & ", maxestep=",dtset%maxestep 2229 2230 if (efield_test_cart(ii) > efield_old_cart(ii)) then 2231 efield_test_cart(ii) = efield_old_cart(ii) + dabs(dtset%maxestep) 2232 else 2233 efield_test_cart(ii) = efield_old_cart(ii) - dabs(dtset%maxestep) 2234 end if 2235 end if 2236 end do 2237 2238 dtset%efield(:) = efield_test_cart(:) 2239 2240 ! !write the field parameters: D, E, P, d, e, p, dbar, ebar, pbar 2241 write(message,'(a,a)') ch10, 'scfcv: Constant reduced ebar and d-field - updating E-field:' 2242 call wrtout(std_out,message,'COLL') 2243 call prtefield(dtset,dtefield,std_out,rprimd) 2244 if(dtset%prtvol>=10)then 2245 call wrtout(ab_out,message,'COLL') 2246 call prtefield(dtset,dtefield,ab_out,rprimd) 2247 end if 2248 2249 2250 ! ! need to update dtset%efield_dot(:) with new value 2251 ! ! This needs to be deleted when efield_dot is deleted 2252 dtefield%efield_dot(1) = dot_product(dtset%efield(:),rprimd(:,1)) 2253 dtefield%efield_dot(2) = dot_product(dtset%efield(:),rprimd(:,2)) 2254 dtefield%efield_dot(3) = dot_product(dtset%efield(:),rprimd(:,3)) 2255 2256 else 2257 2258 write(message,'(a,a)') ch10, 'scfcv: Constant reduced ebar and d-field - Pre E-field:' 2259 call wrtout(std_out,message,'COLL') 2260 call prtefield(dtset,dtefield,std_out,rprimd) 2261 if(dtset%prtvol>=10)then 2262 call wrtout(ab_out,message,'COLL') 2263 call prtefield(dtset,dtefield,ab_out,rprimd) 2264 end if 2265 2266 end if ! scfcv_step > 1 2267 2268 efield_old_cart(:)=dtset%efield(:) 2269 red_efield2_old(:)=red_efield2(:) 2270 2271 end if ! berryopt ==17 2272 2273 end if ! end efield .and. scfcv_level 1 tasks 2274 2275 !deallocate cprj 2276 if ( efield .and. psps%usepaw == 1) then 2277 call pawcprj_free(cprj) 2278 end if 2279 ABI_FREE(cprj) 2280 2281 end subroutine update_e_field_vars