TABLE OF CONTENTS
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-2018 ABINIT group (MVeithen). This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
dtset <type(dataset_type)> = all input variables in this dataset gmet(3,3) = reciprocal space metric tensor in bohr**-2 gprimd(3,3) = primitive translations in recip space kg(3,mpw*mkmem) = reduced (integer) coordinates of G vecs in basis sphere mband = maximum number of bands mkmem = maximum number of k-points in core memory mpw = maximum number of plane waves natom = number of atoms in unit cell nkpt = number of k points npwarr(nkpt) = number of planewaves in basis and boundary at this k point nsppol = 1 for unpolarized, 2 for spin-polarized nsym = number of symmetry operations ntypat = number of types of atoms in unit cell occ(mband*nkpt*nsppol) = occup number for each band at each k point pawang <type(pawang_type)>=paw angular mesh and related data pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data psps <type(pseudopotential_type)>=variables related to pseudopotentials rprimd(3,3) = dimensional primitive vectors symrec(3,3,nsym) = symmetries in reciprocal space in terms of reciprocal space primitive translations typat = typat(natom) list of atom types usepaw = flag for PAW (1 PAW, 0 NCPP) xred(3,natom) = location of atoms in reduced units
OUTPUT
dtefield <type(efield_type)> = variables related to Berry phase calculations pwind(pwind_alloc,2,3) = array used to compute the overlap matrix smat between k-points k and k +- dk where dk is parallel to the direction idir jpw = pwind(ipw,ifor,idir) * ipw = index of plane wave vector G for a given k-point k * ifor = 1: k + dk 2: k - dk * idir = direction of the polarization/ddk calculation [dk(idir) is the only non-zero element of dk(:)] * jpw = index of plane wave vector G (+dG) at k +- dk where dG is a shift of one reciprocal lattice vector (required to close the strings of k-points using the periodic gauge condition) In case a G-vector of the basis sphere of plane waves at k does not belong to the basis sphere of plane waves at k+dk, jpw = 0. pwind_alloc = first dimension of pwind and pwnsfac pwnsfac(2,pwind_alloc) = phase factors for non-symmorphic translations
SIDE EFFECTS
mpi_enreg = informations about MPI parallelization
kptdstrb(nproc,nneighbour,fmkmem_max*nsppol) : Array required
by berryphase_new.f for MPI // over k-points. Defined
for k-points in the fBZ
but for k-points in the iBZ. Used by vtorho.f
nproc = number of cpus
nneighbour = number of neighbours for each k-point (= 6)
PARENTS
init_e_field_vars
CHILDREN
expibi,kpgsph,listkk,metric,pawcprj_alloc,pawcprj_getdim,qijb_kk setsymrhoij,smpbz,symatm,timab,wrtout,xmpi_max,xmpi_sum
SOURCE
82 #if defined HAVE_CONFIG_H 83 #include "config.h" 84 #endif 85 86 #include "abi_common.h" 87 88 subroutine initberry(dtefield,dtset,gmet,gprimd,kg,mband,& 89 & mkmem,mpi_enreg,mpw,natom,nkpt,npwarr,nsppol,& 90 & nsym,ntypat,occ,pawang,pawrad,pawtab,psps,& 91 & pwind,pwind_alloc,pwnsfac,& 92 & rprimd,symrec,typat,usepaw,xred) 93 94 use defs_basis 95 use defs_datatypes 96 use defs_abitypes 97 use m_profiling_abi 98 use m_errors 99 use m_xmpi 100 use m_efield 101 102 use m_geometry,only : metric 103 use m_fftcore, only : kpgsph 104 use m_pawang, only : pawang_type 105 use m_pawrad, only : pawrad_type 106 use m_pawtab, only : pawtab_type 107 use m_pawcprj, only : pawcprj_alloc, pawcprj_getdim 108 109 !This section has been created automatically by the script Abilint (TD). 110 !Do not modify the following lines by hand. 111 #undef ABI_FUNC 112 #define ABI_FUNC 'initberry' 113 use interfaces_14_hidewrite 114 use interfaces_18_timing 115 use interfaces_32_util 116 use interfaces_41_geometry 117 use interfaces_56_recipspace 118 use interfaces_65_paw 119 !End of the abilint section 120 121 implicit none 122 123 !Arguments ------------------------------------ 124 !scalars 125 integer,intent(in) :: mband,mkmem,mpw,natom,nkpt,nsppol,nsym,ntypat,usepaw 126 integer,intent(out) :: pwind_alloc 127 type(MPI_type),intent(inout) :: mpi_enreg 128 type(dataset_type),intent(inout) :: dtset 129 type(efield_type),intent(inout) :: dtefield !vz_i 130 type(pawang_type),intent(in) :: pawang 131 type(pseudopotential_type),intent(in) :: psps 132 !arrays 133 integer,intent(in) :: kg(3,mpw*mkmem),npwarr(nkpt) 134 integer,intent(in) :: symrec(3,3,nsym),typat(natom) 135 integer,pointer :: pwind(:,:,:) 136 real(dp),intent(in) :: gmet(3,3),gprimd(3,3),occ(mband*nkpt*nsppol) 137 real(dp),intent(in) :: rprimd(3,3),xred(3,natom) 138 real(dp),pointer :: pwnsfac(:,:) 139 type(pawrad_type),intent(in) :: pawrad(ntypat) 140 type(pawtab_type),intent(in) :: pawtab(ntypat) 141 142 !Local variables------------------------------- 143 !scalars 144 integer :: exchn2n3d,flag,flag_kpt,fnkpt_computed,iband,icg,icprj 145 integer :: idir,idum,idum1,ierr,ifor,ikg,ikg1,ikpt,ikpt1,ikpt1f 146 integer :: ikpt1i,ikpt2,ikpt_loc,ikptf,ikpti,ikstr,index,ineigh,ipw,ipwnsfac 147 integer :: isppol,istr,istwf_k,isym,isym1,itrs,itypat,iunmark,jpw,klmn,lmax,lmn2_size_max 148 integer :: me,me_g0,mkmem_,my_nspinor,nband_k,mband_occ_k,ncpgr,nkstr,nproc,npw_k,npw_k1,spaceComm 149 integer :: option, brav, mkpt, nkptlatt 150 integer :: jstr,ii,jj,isign 151 integer :: dk_flag, coord1, coord2 152 integer :: mult 153 real(dp) :: c1,ecut_eff,eg,eg_ev,rdum,diffk1,diffk2,diffk3 154 real(dp) :: dist_, max_dist, last_dist, dist,kpt_shifted1,kpt_shifted2,kpt_shifted3 155 real(dp) :: gprimdlc(3,3),rmetllc(3,3),gmetllc(3,3),ucvol_local 156 ! gprimd(3,3) = inverse of rprimd 157 ! rmetlcl(3,3)=real-space metric (same as rmet in metric.F90) 158 ! gmetlcl(3,3)= same as gmet in metric.F90 159 ! ucvol = volume of the unit cell in Bohr**3 160 161 character(len=500) :: message 162 logical :: calc_epaw3_force,calc_epaw3_stress,fieldflag 163 !arrays 164 integer :: dg(3),iadum(3),iadum1(3),neigh(6) 165 integer,allocatable :: dimlmn(:),kg1_k(:,:),kpt_mark(:),nattyp_dum(:) 166 real(dp) :: diffk(3),dk(3),dum33(3,3),eg_dir(3) 167 real(dp) :: kpt1(3) 168 real(dp) :: delta_str3(2), dstr(2),dk_str(2,2,3) 169 real(dp) :: tsec(2) 170 real(dp),allocatable :: calc_expibi(:,:),calc_qijb(:,:,:),spkpt(:,:) 171 172 ! ************************************************************************* 173 174 DBG_ENTER("COLL") 175 176 call timab(1001,1,tsec) 177 call timab(1002,1,tsec) 178 179 !save the current value of berryopt 180 dtefield%berryopt = dtset%berryopt 181 !save the current value of nspinor 182 dtefield%nspinor = dtset%nspinor 183 184 !---------------------------------------------------------------------------- 185 !-------------------- Obtain k-point grid in the full BZ -------------------- 186 !---------------------------------------------------------------------------- 187 188 if(dtset%kptopt==1 .or. dtset%kptopt==2 .or. dtset%kptopt==4)then 189 ! Compute the number of k points in the G-space unit cell 190 nkptlatt=dtset%kptrlatt(1,1)*dtset%kptrlatt(2,2)*dtset%kptrlatt(3,3) & 191 & +dtset%kptrlatt(1,2)*dtset%kptrlatt(2,3)*dtset%kptrlatt(3,1) & 192 & +dtset%kptrlatt(1,3)*dtset%kptrlatt(2,1)*dtset%kptrlatt(3,2) & 193 & -dtset%kptrlatt(1,2)*dtset%kptrlatt(2,1)*dtset%kptrlatt(3,3) & 194 & -dtset%kptrlatt(1,3)*dtset%kptrlatt(2,2)*dtset%kptrlatt(3,1) & 195 & -dtset%kptrlatt(1,1)*dtset%kptrlatt(2,3)*dtset%kptrlatt(3,2) 196 197 ! Call smpbz to obtain the list of k-point in the full BZ - without symmetry reduction 198 option = 0 199 brav = 1 200 mkpt=nkptlatt*dtset%nshiftk 201 ABI_ALLOCATE(spkpt,(3,mkpt)) 202 call smpbz(1,ab_out,dtset%kptrlatt,mkpt,fnkpt_computed,dtset%nshiftk,option,dtset%shiftk,spkpt) 203 dtefield%fnkpt = fnkpt_computed 204 ABI_ALLOCATE(dtefield%fkptns,(3,dtefield%fnkpt)) 205 dtefield%fkptns(:,:)=spkpt(:,1:dtefield%fnkpt) 206 ABI_DEALLOCATE(spkpt) 207 else if(dtset%kptopt==3.or.dtset%kptopt==0)then 208 dtefield%fnkpt=nkpt 209 ABI_ALLOCATE(dtefield%fkptns,(3,dtefield%fnkpt)) 210 dtefield%fkptns(1:3,1:dtefield%fnkpt)=dtset%kpt(1:3,1:dtefield%fnkpt) 211 if(dtset%kptopt==0)then 212 write(message,'(10a)') ch10,& 213 & ' initberry : WARNING -',ch10,& 214 & ' you have defined manually the k-point grid with kptopt = 0',ch10,& 215 & ' the berry phase calculation works only with a regular k-points grid,',ch10,& 216 & ' abinit doesn''t check if your grid is regular...' 217 call wrtout(std_out,message,'PERS') 218 end if 219 end if 220 221 !call listkk to get mapping from FBZ to IBZ 222 rdum=1.0d-5 ! cutoff distance to decide when two k points match 223 ABI_ALLOCATE(dtefield%indkk_f2ibz,(dtefield%fnkpt,6)) 224 225 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 226 227 !ji: The following may need modification in the future 228 !**** no spin-polarization doubling ; allow use of time reversal symmetry **** 229 230 !Here is original call 231 ! 232 !call listkk(rdum,gmet,dtefield%indkk_f2ibz,dtset%kptns,dtefield%fkptns,nkpt,& 233 !& dtefield%fnkpt,dtset%nsym,1,dtset%symafm,dtset%symrel,1) 234 235 call timab(1002,2,tsec) 236 call timab(1003,1,tsec) 237 238 call listkk(rdum,gmet,dtefield%indkk_f2ibz,dtset%kptns,dtefield%fkptns,nkpt,& 239 & dtefield%fnkpt,dtset%nsym,1,dtset%symafm,symrec,1,use_symrec=.True.) 240 241 call timab(1003,2,tsec) 242 call timab(1004,1,tsec) 243 244 !Construct i2fbz and f2ibz 245 ABI_ALLOCATE(dtefield%i2fbz,(nkpt)) 246 idum=0 247 do ikpt=1,dtefield%fnkpt 248 if (dtefield%indkk_f2ibz(ikpt,2)==1 .and. & 249 & dtefield%indkk_f2ibz(ikpt,6) == 0 .and. & 250 & maxval(abs(dtefield%indkk_f2ibz(ikpt,3:5))) == 0 ) then 251 dtefield%i2fbz(dtefield%indkk_f2ibz(ikpt,1))=ikpt 252 idum=idum+1 253 end if 254 end do 255 if (idum/=nkpt)then 256 message = ' Found wrong number of k-points in IBZ' 257 MSG_ERROR(message) 258 end if 259 260 !set flags for fields, forces, stresses 261 fieldflag = ( (dtset%berryopt== 4) .or. (dtset%berryopt== 6) .or. (dtset%berryopt== 7) & 262 & .or. (dtset%berryopt==14) .or. (dtset%berryopt==16) .or. (dtset%berryopt==17) ) 263 ! following two flags activates computation of projector gradient contributions to force and 264 ! stress in finite field PAW calculations 265 calc_epaw3_force = (fieldflag .and. usepaw == 1 .and. dtset%optforces /= 0) 266 calc_epaw3_stress = (fieldflag .and. usepaw == 1 .and. dtset%optstress /= 0) 267 268 269 270 !---------------------------------------------------------------------------- 271 !------------- Allocate PAW space if necessary ------------------------------ 272 !---------------------------------------------------------------------------- 273 274 if (usepaw == 1) then 275 276 dtefield%usepaw = usepaw 277 dtefield%natom = natom 278 dtefield%my_natom = mpi_enreg%my_natom 279 280 ABI_ALLOCATE(dtefield%lmn_size,(ntypat)) 281 ABI_ALLOCATE(dtefield%lmn2_size,(ntypat)) 282 do itypat = 1, ntypat 283 dtefield%lmn_size(itypat) = pawtab(itypat)%lmn_size 284 dtefield%lmn2_size(itypat) = pawtab(itypat)%lmn2_size 285 end do 286 287 lmn2_size_max = psps%lmnmax*(psps%lmnmax+1)/2 288 dtefield%lmn2max = lmn2_size_max 289 290 ! expibi and qijb_kk are NOT parallelized over atoms 291 ! this may change in the future (JZwanziger 18 March 2014) 292 ABI_ALLOCATE(dtefield%qijb_kk,(2,lmn2_size_max,dtefield%natom,3)) 293 ABI_ALLOCATE(dtefield%expibi,(2,dtefield%natom,3)) 294 dtefield%has_expibi = 1 295 dtefield%has_qijb = 1 296 297 if ( fieldflag .and. dtefield%has_rij==0) then 298 lmn2_size_max = psps%lmnmax*(psps%lmnmax+1)/2 299 ABI_ALLOCATE(dtefield%rij,(lmn2_size_max,ntypat,3)) 300 dtefield%has_rij = 1 301 end if 302 303 ! additional F3-type force term for finite electric field with PAW. Same term 304 ! might also apply for other displacement-type field calculations, but not sure yet 305 ! JZwanziger 4 April 2014 306 if ( calc_epaw3_force ) then 307 ABI_ALLOCATE(dtefield%epawf3,(dtefield%natom,3,3)) 308 dtefield%has_epawf3 = 1 309 end if 310 if ( calc_epaw3_stress ) then 311 ABI_ALLOCATE(dtefield%epaws3,(dtefield%natom,3,6)) 312 dtefield%has_epaws3 = 1 313 end if 314 315 ncpgr = 0 316 if ( fieldflag .and. dtefield%usecprj == 0) then 317 ABI_ALLOCATE(dimlmn,(natom)) 318 call pawcprj_getdim(dimlmn,natom,nattyp_dum,ntypat,typat,pawtab,'R') 319 ! allocate space for cprj at kpts in BZ (IBZ or FBZ) 320 ABI_DATATYPE_ALLOCATE(dtefield%cprj,(natom, mband*dtset%nspinor*dtset%nkpt*nsppol)) 321 ! write(std_out,*) "initberry alloc of cprj ", shape(dtefield%cprj) 322 if (calc_epaw3_force .and. .not. calc_epaw3_stress) ncpgr = 3 323 if (.not. calc_epaw3_force .and. calc_epaw3_stress) ncpgr = 6 324 if (calc_epaw3_force .and. calc_epaw3_stress) ncpgr = 9 325 call pawcprj_alloc(dtefield%cprj,ncpgr,dimlmn) 326 dtefield%usecprj = 1 327 ABI_DEALLOCATE(dimlmn) 328 end if 329 330 ABI_ALLOCATE(dtefield%cprjindex,(nkpt,nsppol)) 331 dtefield%cprjindex(:,:) = 0 332 333 if (dtset%kptopt /= 3) then 334 ABI_ALLOCATE(dtefield%atom_indsym,(4,nsym,natom)) 335 call symatm(dtefield%atom_indsym,natom,nsym,symrec,dtset%tnons,tol8,typat,xred) 336 lmax = psps%mpsang - 1 337 ABI_ALLOCATE(dtefield%zarot,(2*lmax+1,2*lmax+1,lmax+1,nsym)) 338 call setsymrhoij(gprimd,lmax,nsym,1,rprimd,symrec,dtefield%zarot) 339 dtefield%nsym = nsym 340 dtefield%lmax = lmax 341 dtefield%lmnmax = psps%lmnmax 342 end if 343 344 end if 345 346 !------------------------------------------------------------------------------ 347 !------------------- Compute variables related to MPI // ---------------------- 348 !------------------------------------------------------------------------------ 349 spaceComm=mpi_enreg%comm_cell 350 nproc=xmpi_comm_size(spaceComm) 351 me=xmpi_comm_rank(spaceComm) 352 353 if (nproc==1) then 354 dtefield%fmkmem = dtefield%fnkpt 355 dtefield%fmkmem_max = dtefield%fnkpt 356 dtefield%mkmem_max = nkpt 357 else 358 dtefield%fmkmem = 0 359 do ikpt = 1, dtefield%fnkpt 360 ikpti = dtefield%indkk_f2ibz(ikpt,1) 361 nband_k = dtset%nband(ikpti) 362 if (.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpti,1,nband_k,-1,me))) & 363 & dtefield%fmkmem = dtefield%fmkmem + 1 364 end do 365 ! Maximum value of mkmem and fmkmem 366 call xmpi_max(dtefield%fmkmem,dtefield%fmkmem_max,spaceComm,ierr) 367 ! I have to use the dummy variable mkmem_ because 368 ! mkmem is declared as intent(in) while the first 369 ! argument of xmpi_max must be intent(inout) 370 mkmem_ = mkmem 371 call xmpi_max(mkmem_,dtefield%mkmem_max,spaceComm,ierr) 372 end if 373 374 ABI_ALLOCATE(mpi_enreg%kpt_loc2fbz_sp,(0:nproc-1,1:dtefield%fmkmem_max*nsppol, 1:2)) 375 ABI_ALLOCATE(mpi_enreg%kpt_loc2ibz_sp,(0:nproc-1,1:dtefield%mkmem_max*nsppol, 1:2)) 376 ABI_ALLOCATE(mpi_enreg%kptdstrb,(nproc,6,dtefield%fmkmem_max*nsppol*2)) 377 ABI_ALLOCATE(mpi_enreg%mkmem,(0:nproc-1)) 378 mpi_enreg%kpt_loc2fbz_sp(:,:,:) = 0 379 mpi_enreg%kpt_loc2ibz_sp(:,:,:) = 0 380 mpi_enreg%kptdstrb(:,:,:) = 0 381 mpi_enreg%mkmem(:) = 0 382 383 if (fieldflag) then 384 ABI_ALLOCATE(dtefield%cgqindex,(3,6,nkpt*nsppol)) 385 ABI_ALLOCATE(dtefield%nneigh,(nkpt)) 386 dtefield%cgqindex(:,:,:) = 0 ; dtefield%nneigh(:) = 0 387 end if 388 389 pwind_alloc = mpw*dtefield%fmkmem_max 390 ABI_ALLOCATE(pwind,(pwind_alloc,2,3)) 391 ABI_ALLOCATE(pwnsfac,(2,pwind_alloc)) 392 393 !------------------------------------------------------------------------------ 394 !---------------------- Compute efield_type variables ------------------------- 395 !------------------------------------------------------------------------------ 396 397 !Initialization of efield_type variables 398 mult=dtset%useria+1 399 dtefield%efield_dot(:) = zero 400 dtefield%dkvecs(:,:) = zero 401 dtefield%maxnstr = 0 ; dtefield%maxnkstr = 0 402 dtefield%nstr(:) = 0 ; dtefield%nkstr(:) = 0 403 ABI_ALLOCATE(dtefield%ikpt_dk,(dtefield%fnkpt,2,3)) 404 ABI_ALLOCATE(dtefield%cgindex,(nkpt,nsppol)) 405 ABI_ALLOCATE(dtefield%kgindex,(nkpt)) 406 ABI_ALLOCATE(dtefield%fkgindex,(dtefield%fnkpt)) 407 dtefield%ikpt_dk(:,:,:) = 0 408 dtefield%cgindex(:,:) = 0 409 dtefield%mband_occ = 0 410 ABI_ALLOCATE(dtefield%nband_occ,(nsppol)) 411 dtefield%kgindex(:) = 0 412 dtefield%fkgindex(:) = 0 413 414 if (fieldflag) then 415 dtset%rfdir(1:3) = 1 416 end if 417 418 419 !Compute spin degeneracy 420 if (nsppol == 1 .and. dtset%nspinor == 1) then 421 dtefield%sdeg = two 422 else if (nsppol == 2 .or. my_nspinor == 2) then 423 dtefield%sdeg = one 424 end if 425 426 !Compute the number of occupied bands and check that 427 !it is the same for each k-point 428 429 index = 0 430 do isppol = 1, nsppol 431 dtefield%nband_occ(isppol) = 0 432 do ikpt = 1, nkpt 433 434 mband_occ_k = 0 435 nband_k = dtset%nband(ikpt + (isppol - 1)*nkpt) 436 437 do iband = 1, nband_k 438 index = index + 1 439 if (abs(occ(index) - dtefield%sdeg) < tol8) mband_occ_k = mband_occ_k + 1 440 end do 441 442 if (fieldflag) then 443 if (nband_k /= mband_occ_k) then 444 write(message,'(a,a,a)')& 445 & ' In a finite electric field, nband must be equal ',ch10,& 446 & ' to the number of valence bands.' 447 MSG_ERROR(message) 448 end if 449 end if 450 451 if (ikpt > 1) then 452 if (dtefield%nband_occ(isppol) /= mband_occ_k) then 453 message = "The number of valence bands is not the same for every k-point of present spin channel" 454 MSG_ERROR(message) 455 end if 456 else 457 dtefield%mband_occ = max(dtefield%mband_occ, mband_occ_k) 458 dtefield%nband_occ(isppol) = mband_occ_k 459 end if 460 461 end do ! close loop over ikpt 462 end do ! close loop over isppol 463 464 if (fieldflag) then 465 ABI_ALLOCATE(dtefield%smat,(2,dtefield%mband_occ,dtefield%mband_occ,nkpt*nsppol,2,3)) 466 467 dtefield%smat(:,:,:,:,:,:) = zero 468 end if 469 470 ABI_ALLOCATE(dtefield%sflag,(dtefield%mband_occ,nkpt*nsppol,2,3)) 471 dtefield%sflag(:,:,:,:) = 0 472 473 !Compute the location of each wavefunction 474 475 icg = 0 476 icprj = 0 477 !ikg = 0 478 do isppol = 1, nsppol 479 do ikpt = 1, nkpt 480 481 nband_k = dtset%nband(ikpt + (isppol-1)*nkpt) 482 483 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) cycle 484 485 dtefield%cgindex(ikpt,isppol) = icg 486 npw_k = npwarr(ikpt) 487 icg = icg + npw_k*dtefield%nspinor*nband_k 488 489 if (usepaw == 1) then 490 dtefield%cprjindex(ikpt,isppol) = icprj 491 icprj = icprj + dtefield%nspinor*nband_k 492 end if 493 494 end do 495 end do 496 497 ikg = 0 498 do ikpt = 1, nkpt 499 if ((proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,1,me)).and.& 500 & (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,nsppol,me))) cycle 501 502 npw_k = npwarr(ikpt) 503 dtefield%kgindex(ikpt) = ikg 504 ikg = ikg + npw_k 505 end do 506 507 !Need to use dtset%red_efieldbar in the whole code 508 !Compute the reciprocal lattice coordinates of the electric field 509 if (fieldflag) then 510 511 call metric(gmetllc,gprimdlc,-1,rmetllc,rprimd,ucvol_local) 512 513 if (dtset%berryopt == 4 .or. dtset%berryopt == 6 .or. dtset%berryopt == 7) then 514 515 do ii=1,3 516 dtset%red_efieldbar(ii) = dot_product(dtset%efield(:),rprimd(:,ii)) 517 dtefield%efield_dot(ii) = dtset%red_efieldbar(ii) 518 end do 519 520 ! dtefield%efield_dot(1) = dot_product(dtset%efield(:),rprimd(:,1)) 521 ! dtefield%efield_dot(2) = dot_product(dtset%efield(:),rprimd(:,2)) 522 ! dtefield%efield_dot(3) = dot_product(dtset%efield(:),rprimd(:,3)) 523 524 write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,& 525 & ' initberry: Reduced electric field (ebar)',ch10,& 526 & ' red_efieldbar(1:3) = ',dtset%red_efieldbar(1:3),ch10 527 call wrtout(std_out,message,'COLL') 528 529 end if 530 531 if (dtset%berryopt == 6 .or. dtset%berryopt ==7 ) then 532 533 do ii=1,3 534 dtset%red_dfield(ii)= (dot_product(dtset%dfield(:),gprimdlc(:,ii)))*ucvol_local/(4.d0*pi) 535 end do 536 537 write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,& 538 & ' initberry: Reduced electric displacement field',ch10,& 539 & ' red_dfield(1:3) = ',dtset%red_dfield(1:3),ch10 540 call wrtout(std_out,message,'COLL') 541 542 end if 543 544 545 if ( dtset%berryopt == 14 ) then 546 ! transfer to unreduced electric field. 547 do idir=1,3 548 dtset%efield(idir)= dot_product(dtset%red_efieldbar(:),gprimdlc(:,idir)) 549 dtefield%efield_dot(idir) = dtset%red_efieldbar(idir) 550 ! dtefield%efield2(idir)=dtset%red_efieldbar(idir) 551 end do 552 553 ! dtefield%efield_dot(1) = dtset%red_efieldbar(1) 554 ! dtefield%efield_dot(2) = dtset%red_efieldbar(2) 555 ! dtefield%efield_dot(3) = dtset%red_efieldbar(3) 556 557 write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,& 558 & ' initberry: Unreduced electric field (a.u.)',ch10,& 559 & ' efield(1:3) = ',dtset%efield(1:3),ch10 560 call wrtout(std_out,message,'COLL') 561 562 write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,& 563 & ' initberry: Reduced electric field (ebar)',ch10,& 564 & ' red_efieldbar(1:3) = ',dtset%red_efieldbar(1:3),ch10 565 call wrtout(std_out,message,'COLL') 566 567 end if 568 569 570 if ( dtset%berryopt == 16 ) then 571 572 ! to calculate D 573 do ii=1,3 574 dtset%dfield(ii) =(4*pi/ucvol_local)*dot_product(dtset%red_dfield(:),rprimd(:,ii)) 575 end do 576 577 do idir=1,3 578 dtset%efield(idir)= (4*pi/ucvol_local)*dot_product(dtset%red_efield(:),rprimd(:,idir)) 579 end do 580 581 do idir=1,3 582 dtset%red_efieldbar(idir)= (4*pi/ucvol_local)*dot_product(dtset%red_efield(:),rmetllc(:,idir)) 583 dtefield%efield_dot(idir) = dtset%red_efieldbar(idir) 584 end do 585 586 ! dtefield%efield_dot(1) = dtset%red_efieldbar(1) 587 ! dtefield%efield_dot(2) = dtset%red_efieldbar(2) 588 ! dtefield%efield_dot(3) = dtset%red_efieldbar(3) 589 590 591 write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,& 592 & ' initberry: Unreduced electric displacement field (a.u.)',ch10,& 593 & ' dfield(1:3) = ',dtset%dfield(1:3),ch10 594 call wrtout(std_out,message,'COLL') 595 596 write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,& 597 & ' initberry: Unreduced electric field (a.u.)',ch10,& 598 & ' efield(1:3) = ',dtset%efield(1:3),ch10 599 call wrtout(std_out,message,'COLL') 600 601 write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,& 602 & ' initberry: Reduced electric field (ebar)',ch10,& 603 & ' red_efieldbar(1:3) = ',dtset%red_efieldbar(1:3),ch10 604 call wrtout(std_out,message,'COLL') 605 606 end if 607 608 if ( dtset%berryopt ==17) then 609 610 ! to calculate D 611 612 do idir=1,3 613 dtset%efield(idir)= dot_product(dtset%red_efieldbar(:),gprimdlc(:,idir)) ! from ebar 614 dtset%dfield(idir) =(4*pi/ucvol_local)*dot_product(dtset%red_dfield(:),rprimd(:,idir)) 615 ! dtset%red_efield(idir) = (ucvol_local/(4*pi))*dot_product(dtset%red_efieldbar(:),gmetllc(:,idir)) 616 dtefield%efield_dot(idir) = dtset%red_efieldbar(idir) 617 end do 618 619 write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,& 620 & ' initberry: Reduced electric field (ebar)',ch10,& 621 & ' red_efieldbar(1:3) = ',dtset%red_efieldbar(1:3),ch10 622 call wrtout(std_out,message,'COLL') 623 624 625 write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,& 626 & ' initberry: Unreduced electric field (a.u.)',ch10,& 627 & ' efield(1:3) = ',dtset%efield(1:3),ch10 628 call wrtout(std_out,message,'COLL') 629 630 write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,& 631 & ' initberry: Reduced electric displacement field (a.u.)',ch10,& 632 & ' red_dfield(1:3) = ',dtset%red_dfield(1:3),ch10 633 call wrtout(std_out,message,'COLL') 634 635 write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,& 636 & ' initberry: Unreduced electric displacement field (a.u.)',ch10,& 637 & ' dfield(1:3) = ',dtset%dfield(1:3),ch10 638 call wrtout(std_out,message,'COLL') 639 640 641 end if 642 643 644 645 end if 646 647 call timab(1004,2,tsec) 648 649 !------------------------------------------------------------------------------ 650 !---------------------- Compute dk -------------------------------------------- 651 !------------------------------------------------------------------------------ 652 653 call timab(1005,1,tsec) 654 655 do idir = 1, 3 656 657 if (dtset%rfdir(idir) == 1) then 658 659 ! Compute dk(:), the vector between a k-point and its nearest 660 ! neighbour along the direction idir 661 662 dk(:) = zero 663 dk(idir) = 1._dp ! 1 mean there is no other k-point un the direction idir 664 do ikpt = 2, dtefield%fnkpt 665 diffk(:) = abs(dtefield%fkptns(:,ikpt) - dtefield%fkptns(:,1)) 666 if ((diffk(1) < dk(1)+tol8).and.(diffk(2) < dk(2)+tol8).and.& 667 & (diffk(3) < dk(3)+tol8)) dk(:) = diffk(:) 668 end do 669 dtefield%dkvecs(:,idir) = dk(:) 670 ! DEBUG 671 ! write(std_out,*)' initberry : idir, dk', idir, dk 672 ! ENDDEBUG 673 674 ! For each k point, find k_prim such that k_prim= k + dk mod(G) 675 ! where G is a vector of the reciprocal lattice 676 677 do ikpt = 1, dtefield%fnkpt 678 679 ! First k+dk, then k-dk 680 do isign=-1,1,2 681 kpt_shifted1=dtefield%fkptns(1,ikpt)- isign*dk(1) 682 kpt_shifted2=dtefield%fkptns(2,ikpt)- isign*dk(2) 683 kpt_shifted3=dtefield%fkptns(3,ikpt)- isign*dk(3) 684 ! Note that this is still a order fnkpt**2 algorithm. 685 ! It is possible to implement a order fnkpt algorithm, see listkk.F90. 686 do ikpt1 = 1, dtefield%fnkpt 687 diffk1=dtefield%fkptns(1,ikpt1) - kpt_shifted1 688 if(abs(diffk1-nint(diffk1))>tol8)cycle 689 diffk2=dtefield%fkptns(2,ikpt1) - kpt_shifted2 690 if(abs(diffk2-nint(diffk2))>tol8)cycle 691 diffk3=dtefield%fkptns(3,ikpt1) - kpt_shifted3 692 if(abs(diffk3-nint(diffk3))>tol8)cycle 693 dtefield%ikpt_dk(ikpt,(isign+3)/2,idir) = ikpt1 694 exit 695 end do ! ikpt1 696 end do ! isign 697 698 ! OLD CODING 699 ! First: k + dk 700 ! do ikpt1 = 1, dtefield%fnkpt 701 ! diffk(:) = abs(dtefield%fkptns(:,ikpt1) - & 702 ! & dtefield%fkptns(:,ikpt) - dk(:)) 703 ! if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then 704 ! dtefield%ikpt_dk(ikpt,1,idir) = ikpt1 705 ! exit 706 ! end if 707 ! end do 708 709 ! Second: k - dk 710 ! do ikpt1 = 1, dtefield%fnkpt 711 ! diffk(:) = abs(dtefield%fkptns(:,ikpt1) - & 712 ! & dtefield%fkptns(:,ikpt) + dk(:)) 713 ! if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then 714 ! dtefield%ikpt_dk(ikpt,2,idir) = ikpt1 715 ! exit 716 ! end if 717 ! end do 718 719 end do ! ikpt 720 721 ! Find the string length, starting from k point 1 722 ! (all strings must have the same number of points) 723 724 nkstr = 1 725 ikpt1 = 1 726 do ikpt = 1, dtefield%fnkpt 727 ikpt1 = dtefield%ikpt_dk(ikpt1,1,idir) 728 if (ikpt1 == 1) exit 729 nkstr = nkstr + 1 730 end do 731 732 ! Check that the string length is a divisor of nkpt 733 if(mod(dtefield%fnkpt,nkstr) /= 0) then 734 write(message,'(a,i5,a,i7)')& 735 & ' The string length = ',nkstr,& 736 & ', is not a divisor of fnkpt =',dtefield%fnkpt 737 MSG_BUG(message) 738 end if 739 740 dtefield%nkstr(idir) = nkstr 741 dtefield%nstr(idir) = dtefield%fnkpt/nkstr 742 743 end if ! dtset%rfdir(idir) == 1 744 745 write(message,'(a,i1,a,i3,a,i6)')& 746 & ' initberry: for direction ',idir,', nkstr = ',dtefield%nkstr(idir),& 747 & ', nstr = ',dtefield%nstr(idir) 748 call wrtout(std_out,message,'COLL') 749 call wrtout(ab_out,message,'COLL') 750 751 end do ! close loop over idir 752 753 call timab(1005,2,tsec) 754 call timab(1006,1,tsec) 755 756 dtefield%maxnstr = maxval(dtefield%nstr(:)) 757 dtefield%maxnkstr = maxval(dtefield%nkstr(:)) 758 ABI_ALLOCATE(dtefield%idxkstr,(dtefield%maxnkstr,dtefield%maxnstr,3)) 759 dtefield%idxkstr(:,:,:) = 0 760 761 !for the geometry of the string space : 762 ABI_ALLOCATE(dtefield%coord_str,(2,dtefield%maxnstr,3)) 763 ABI_ALLOCATE(dtefield%str_neigh,(-2:2,dtefield%maxnstr,3)) 764 ABI_ALLOCATE(dtefield%strg_neigh,(-2:2,dtefield%maxnstr,2,3)) 765 dtefield%coord_str(:,:,:) = 0.d0 766 dtefield%str_neigh(:,:,:)=0 767 dtefield%strg_neigh(:,:,:,:)=0 768 dtefield%gmet_str(:,:,:)=0.d0 769 770 !------------------------------------------------------------------------------ 771 !---------------------- Build the strings ------------------------------------- 772 !------------------------------------------------------------------------------ 773 774 ABI_ALLOCATE(kpt_mark,(dtefield%fnkpt)) 775 do idir = 1, 3 776 777 if (dtset%rfdir(idir) == 1) then 778 779 iunmark = 1 780 kpt_mark(:) = 0 781 do istr = 1, dtefield%nstr(idir) 782 783 do while(kpt_mark(iunmark) /= 0) 784 iunmark = iunmark + 1 785 end do 786 dtefield%idxkstr(1,istr,idir) = iunmark 787 kpt_mark(iunmark) = 1 788 do ikstr = 2, dtefield%nkstr(idir) 789 ikpt1 = dtefield%idxkstr(ikstr-1,istr,idir) 790 ikpt2 = dtefield%ikpt_dk(ikpt1,1,idir) 791 dtefield%idxkstr(ikstr,istr,idir) = ikpt2 792 kpt_mark(ikpt2) = 1 793 end do 794 795 end do ! istr 796 797 ! compute distance between strings 798 ! compute the metric matrix of the strings space in the direction idir 799 do ii = 1,3 800 do jj = 1,3 801 if (ii<idir.and.jj<idir) dtefield%gmet_str(ii ,jj ,idir) = & 802 & gmet(ii,jj) - gmet(ii,idir)*gmet(jj,idir)/gmet(idir,idir) 803 if (ii<idir.and.jj>idir) dtefield%gmet_str(ii ,jj-1,idir) = & 804 & gmet(ii,jj) - gmet(ii,idir)*gmet(jj,idir)/gmet(idir,idir) 805 if (ii>idir.and.jj<idir) dtefield%gmet_str(ii-1,jj ,idir) = & 806 & gmet(ii,jj) - gmet(ii,idir)*gmet(jj,idir)/gmet(idir,idir) 807 if (ii>idir.and.jj>idir) dtefield%gmet_str(ii-1,jj-1,idir) = & 808 & gmet(ii,jj) - gmet(ii,idir)*gmet(jj,idir)/gmet(idir,idir) 809 end do 810 end do 811 ! DEBUG 812 ! write(std_out,*)'gmet' 813 ! do ii=1,3 814 ! write(std_out,*)gmet(ii,:) 815 ! end do 816 ! write(std_out,*)'gmet_str' 817 ! do ii=1,2 818 ! write(std_out,*)dtefield%gmet_str(ii,:,idir) 819 ! end do 820 ! ENDDEBUG 821 do istr = 1, dtefield%nstr(idir) 822 do ii = 1,3 823 if (ii<idir) dtefield%coord_str(ii,istr,idir)=dtefield%fkptns(ii,dtefield%idxkstr(1,istr,idir)) 824 if (ii>idir) dtefield%coord_str(ii-1,istr,idir)=dtefield%fkptns(ii,dtefield%idxkstr(1,istr,idir)) 825 end do 826 end do 827 828 ! the following is very similar to getshell 829 dist_ = 0._dp 830 do ii = 1,2 831 dist_ = dist_ + dtefield%gmet_str(ii,ii,idir) 832 end do 833 max_dist = 2._dp * dist_ * 2._dp 834 835 dk_str(:,:,idir) = 0._dp 836 last_dist = 0._dp 837 ! ishell = 0 838 ! dtefield%str_neigh(:,:,:) = 0 839 dk_flag = 0 840 do while (dk_flag /= 2) 841 ! Advance shell counter 842 ! ishell = ishell + 1 843 844 ! Search the smallest distance between two strings 845 dist = max_dist 846 do istr = 1,dtefield%nstr(idir) 847 delta_str3(:) = dtefield%coord_str(:,1,idir) - dtefield%coord_str(:,istr,idir) 848 do coord1 = -1,1 !two loop to search also on the border of the BZ 849 do coord2 = -1,1 850 dist_ = 0._dp 851 dstr(:) = delta_str3(:) - nint(delta_str3(:)) 852 dstr(1) = dstr(1) + real(coord1,dp) 853 dstr(2) = dstr(2) + real(coord2,dp) 854 do ii = 1,2 855 do jj = 1,2 856 dist_ = dist_ + dstr(ii)*dtefield%gmet_str(ii,jj,idir)*dstr(jj) 857 end do 858 end do 859 if ((dist_ < dist).and.(dist_ - last_dist > tol8)) then 860 dist = dist_ 861 end if 862 end do 863 end do 864 end do 865 866 last_dist = dist 867 868 ! search the connecting vectors for that distance 869 do istr = 1,dtefield%nstr(idir) 870 delta_str3(:) = dtefield%coord_str(:,istr,idir) - dtefield%coord_str(:,1,idir) 871 do coord1 = -1,1 872 do coord2 = -1,1 873 dist_ = 0._dp 874 dstr(:) = delta_str3(:) - nint(delta_str3(:)) 875 dstr(1) = dstr(1) + real(coord1,dp) 876 dstr(2) = dstr(2) + real(coord2,dp) 877 do ii = 1,2 878 do jj = 1,2 879 dist_ = dist_ + dstr(ii)*dtefield%gmet_str(ii,jj,idir)*dstr(jj) 880 end do 881 end do 882 if (abs(dist_ - dist) < tol8) then 883 if (dk_flag == 0) then 884 dk_str(:,1,idir) = dstr(:) 885 dk_flag = 1 886 ! DEBUG 887 ! write(std_out,'(a,i4,2e15.4)')'1st connect', istr, dstr 888 ! ENDDEBUG 889 elseif (dk_str(1,1,idir)*dstr(2)-dk_str(2,1,idir)*dstr(1) > tol8) then 890 dk_str(:,2,idir) = dstr(:) 891 dk_flag = 2 892 ! DEBUG 893 ! write(std_out,'(a,i4,2e15.4)')'2nd connect', istr, dstr 894 ! ENDDEBUG 895 exit 896 end if 897 end if 898 end do 899 if (dk_flag == 2) exit 900 end do 901 if (dk_flag == 2) exit 902 end do 903 904 end do ! do while 905 906 ! search the two neighbours for each string 907 do istr = 1,dtefield%nstr(idir) 908 dtefield%str_neigh(0,istr,idir) = istr 909 dtefield%strg_neigh(0,istr,:,idir) = 0 910 do jstr = 1,dtefield%nstr(idir) 911 delta_str3(:) = dtefield%coord_str(:,jstr,idir) - dtefield%coord_str(:,istr,idir) 912 do coord1 = -1,1 913 do coord2 = -1,1 914 dist_ = 0._dp 915 dstr(:) = delta_str3(:) - nint(delta_str3(:)) 916 dstr(1) = dstr(1) + real(coord1,dp) 917 dstr(2) = dstr(2) + real(coord2,dp) 918 do ii = 1,2 919 if (sum(abs(dstr(:)-dk_str(:,ii,idir)))<tol8) then 920 dtefield%str_neigh(ii,istr,idir) = jstr 921 dtefield%strg_neigh(ii,istr,1,idir) = coord1 922 dtefield%strg_neigh(ii,istr,2,idir) = coord2 923 elseif (sum(abs(dstr(:)+dk_str(:,ii,idir)))<tol8) then 924 dtefield%str_neigh(-ii,istr,idir) = jstr 925 dtefield%strg_neigh(-ii,istr,1,idir) = coord1 926 dtefield%strg_neigh(-ii,istr,2,idir) = coord2 927 end if 928 end do 929 end do 930 end do 931 end do 932 end do 933 934 ! DEBUG 935 ! 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) 936 ! write(std_out,*)'istr, neigh1, strg(1,:), neigh2, strg(2,:),neigh-1, strg(-1,:), neigh-2, strg(-2,:)' 937 ! do istr=1,dtefield%nstr(idir) 938 ! write(std_out,'(13i4)')istr, & 939 ! & dtefield%str_neigh(1,istr,idir), dtefield%strg_neigh(1,istr,:,idir),& 940 ! & dtefield%str_neigh(2,istr,idir), dtefield%strg_neigh(2,istr,:,idir),& 941 ! & dtefield%str_neigh(-1,istr,idir), dtefield%strg_neigh(-1,istr,:,idir),& 942 ! & dtefield%str_neigh(-2,istr,idir), dtefield%strg_neigh(-2,istr,:,idir) 943 ! end do 944 ! ENDDEBUG 945 946 947 end if ! rfdir(idir) == 1 948 949 end do ! close loop over idir 950 951 ABI_DEALLOCATE(kpt_mark) 952 953 call timab(1006,2,tsec) 954 call timab(1007,1,tsec) 955 956 !------------------------------------------------------------------------------ 957 !------------ Compute PAW on-site terms if necessary -------------------------- 958 !------------------------------------------------------------------------------ 959 960 if (usepaw == 1 .and. dtefield%has_expibi == 1) then 961 ABI_ALLOCATE(calc_expibi,(2,natom)) 962 do idir = 1, 3 963 dk = dtefield%dkvecs(1:3,idir) 964 calc_expibi = zero 965 call expibi(calc_expibi,dk,natom,xred) 966 dtefield%expibi(1:2,1:natom,idir) = calc_expibi 967 end do 968 ! call expibi(dtefield%expibi,dtefield%dkvecs,natom,xred) 969 dtefield%has_expibi = 2 970 ABI_DEALLOCATE(calc_expibi) 971 end if 972 973 if (usepaw == 1 .and. dtefield%has_qijb == 1) then 974 ABI_ALLOCATE(calc_qijb,(2,dtefield%lmn2max,natom)) 975 976 do idir = 1, 3 977 dk = dtefield%dkvecs(1:3,idir) 978 calc_qijb = zero 979 call qijb_kk(calc_qijb,dk,dtefield%expibi(1:2,1:natom,idir),& 980 & gprimd,dtefield%lmn2max,natom,ntypat,pawang,pawrad,pawtab,typat) 981 dtefield%qijb_kk(1:2,1:dtefield%lmn2max,1:natom,idir) = calc_qijb 982 ! call qijb_kk(dtefield%qijb_kk,dtefield%dkvecs,dtefield%expibi,& 983 ! & gprimd,dtefield%lmn2max,natom,ntypat,pawang,pawrad,pawtab,typat) 984 end do 985 dtefield%has_qijb = 2 986 ABI_DEALLOCATE(calc_qijb) 987 end if 988 989 if (usepaw == 1 .and. dtefield%has_rij == 1) then 990 c1=sqrt(four_pi/three) 991 do itypat = 1, ntypat 992 do klmn = 1, pawtab(itypat)%lmn2_size 993 dtefield%rij(klmn,itypat,1) = c1*pawtab(itypat)%qijl(4,klmn) ! S_{1,1} ~ x 994 dtefield%rij(klmn,itypat,2) = c1*pawtab(itypat)%qijl(2,klmn) ! S_{1,-1} ~ y 995 dtefield%rij(klmn,itypat,3) = c1*pawtab(itypat)%qijl(3,klmn) ! S_{1,0} ~ z 996 end do ! end loop over klmn 997 end do ! end loop over itypat 998 dtefield%has_rij = 2 999 end if ! 1000 1001 call timab(1007,2,tsec) 1002 call timab(1008,1,tsec) 1003 1004 !------------------------------------------------------------------------------ 1005 !------------ Build the array pwind that is needed to compute the ------------- 1006 !------------ overlap matrices at k +- dk ------------- 1007 !------------------------------------------------------------------------------ 1008 1009 ecut_eff = dtset%ecut*(dtset%dilatmx)**2 1010 exchn2n3d = 0 ; istwf_k = 1 ; ikg1 = 0 1011 pwind(:,:,:) = 0 1012 pwnsfac(1,:) = 1.0_dp 1013 pwnsfac(2,:) = 0.0_dp 1014 ABI_ALLOCATE(kg1_k,(3,mpw)) 1015 1016 ipwnsfac = 0 1017 1018 do idir = 1, 3 1019 1020 if (dtset%rfdir(idir) == 1) then 1021 1022 dk(:) = dtefield%dkvecs(:,idir) 1023 1024 do ifor = 1, 2 1025 1026 if (ifor == 2) dk(:) = -1._dp*dk(:) 1027 1028 ! Build pwind and kgindex 1029 ! NOTE: The array kgindex is important for parallel execution. 1030 ! In case nsppol = 2, it may happent that a particular processor 1031 ! treats k-points at different spin polarizations. 1032 ! In this case, it is not possible to address the elements of 1033 ! pwind correctly without making use of the kgindex array. 1034 1035 ikg = 0 ; ikpt_loc = 0 ; isppol = 1 1036 do ikpt = 1, dtefield%fnkpt 1037 1038 ikpti = dtefield%indkk_f2ibz(ikpt,1) 1039 nband_k = dtset%nband(ikpti) 1040 ikpt1f = dtefield%ikpt_dk(ikpt,ifor,idir) 1041 ikpt1i = dtefield%indkk_f2ibz(ikpt1f,1) 1042 1043 if ((proc_distrb_cycle(mpi_enreg%proc_distrb,ikpti,1,nband_k,1,me)).and.& 1044 & (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpti,1,nband_k,nsppol,me))) cycle 1045 1046 ikpt_loc = ikpt_loc + 1 1047 1048 ! Build basis sphere of plane waves for the nearest neighbour of 1049 ! the k-point (important for MPI //) 1050 1051 kg1_k(:,:) = 0 1052 kpt1(:) = dtset%kptns(:,ikpt1i) 1053 call kpgsph(ecut_eff,exchn2n3d,gmet,ikg1,ikpt,istwf_k,kg1_k,kpt1,& 1054 & 1,mpi_enreg,mpw,npw_k1) 1055 me_g0=mpi_enreg%me_g0 1056 1057 1058 ! ji: fkgindex is defined here ! 1059 dtefield%fkgindex(ikpt) = ikg 1060 1061 ! 1062 ! Deal with symmetry transformations 1063 ! 1064 1065 ! bra k-point k(b) and IBZ k-point kIBZ(b) related by 1066 ! k(b) = alpha(b) S(b)^t kIBZ(b) + G(b) 1067 ! where alpha(b), S(b) and G(b) are given by indkk_f2ibz 1068 ! 1069 ! For the ket k-point: 1070 ! k(k) = alpha(k) S(k)^t kIBZ(k) + G(k) - GBZ(k) 1071 ! where GBZ(k) takes k(k) to the BZ 1072 ! 1073 1074 isym = dtefield%indkk_f2ibz(ikpt,2) 1075 isym1 = dtefield%indkk_f2ibz(ikpt1f,2) 1076 1077 ! Construct transformed G vector that enters the matching condition: 1078 ! alpha(k) S(k)^{t,-1} ( -G(b) - GBZ(k) + G(k) ) 1079 1080 dg(:) = -dtefield%indkk_f2ibz(ikpt,3:5) & 1081 & -nint(-dtefield%fkptns(:,ikpt) - dk(:) - tol10 + & 1082 & dtefield%fkptns(:,ikpt1f)) & 1083 & +dtefield%indkk_f2ibz(ikpt1f,3:5) 1084 1085 ! old code 1086 ! iadum(:)=0 1087 ! do idum=1,3 1088 ! iadum(:)=iadum(:)+ symrec(:,idum,isym1)*dg(idum) 1089 ! end do 1090 1091 ! new code 1092 iadum(:) = MATMUL(TRANSPOSE(dtset%symrel(:,:,isym1)),dg(:)) 1093 1094 dg(:) = iadum(:) 1095 1096 if ( dtefield%indkk_f2ibz(ikpt1f,6) == 1 ) dg(:) = -dg(:) 1097 1098 ! Construct S(k)^{t,-1} S(b)^{t} 1099 1100 dum33(:,:) = MATMUL(TRANSPOSE(dtset%symrel(:,:,isym1)),symrec(:,:,isym)) 1101 1102 ! Construct alpha(k) alpha(b) 1103 1104 if (dtefield%indkk_f2ibz(ikpt,6) == dtefield%indkk_f2ibz(ikpt1f,6)) then 1105 itrs=0 1106 else 1107 itrs=1 1108 end if 1109 1110 1111 npw_k = npwarr(ikpti) 1112 ! npw_k1 = npwarr(ikpt1i) 1113 1114 ! loop over bra G vectors 1115 do ipw = 1, npw_k 1116 1117 ! NOTE: the bra G vector is taken for the sym-related IBZ k point, 1118 ! not for the FBZ k point 1119 iadum(:) = kg(:,dtefield%kgindex(ikpti) + ipw) 1120 1121 ! Store non-symmorphic operation phase factor exp[i2\pi \alpha G \cdot t] 1122 1123 if ( ipwnsfac == 0 ) then 1124 ! old code 1125 rdum=0.0_dp 1126 do idum=1,3 1127 rdum=rdum+dble(iadum(idum))*dtset%tnons(idum,isym) 1128 end do 1129 rdum=two_pi*rdum 1130 if ( dtefield%indkk_f2ibz(ikpt,6) == 1 ) rdum=-rdum 1131 pwnsfac(1,ikg+ipw) = cos(rdum) 1132 pwnsfac(2,ikg+ipw) = sin(rdum) 1133 ! 1134 ! new code 1135 ! rdum = DOT_PRODUCT(dble(iadum(:)),dtset%tnons(:,isym)) 1136 ! rdum= two_pi*rdum 1137 ! if ( dtefield%indkk_f2ibz(ikpt,6) == 1 ) rdum=-rdum 1138 ! pwnsfac(1,ikg+ipw) = cos(rdum) 1139 ! pwnsfac(2,ikg+ipw) = sin(rdum) 1140 1141 end if 1142 1143 ! to determine r.l.v. matchings, we transformed the bra vector 1144 ! Rotation 1145 iadum1(:)=0 1146 do idum1=1,3 1147 iadum1(:)=iadum1(:)+dum33(:,idum1)*iadum(idum1) 1148 end do 1149 iadum(:)=iadum1(:) 1150 ! Time reversal 1151 if (itrs==1) iadum(:)=-iadum(:) 1152 ! Translation 1153 iadum(:) = iadum(:) + dg(:) 1154 1155 do jpw = 1, npw_k1 1156 iadum1(1:3) = kg1_k(1:3,jpw) 1157 if ( (iadum(1) == iadum1(1)).and. & 1158 & (iadum(2) == iadum1(2)).and. & 1159 & (iadum(3) == iadum1(3)) ) then 1160 pwind(ikg + ipw,ifor,idir) = jpw 1161 ! write(std_out,'(a,2x,3i4,2x,i4)') 'Found !:',iadum1(:),jpw 1162 exit 1163 end if 1164 end do 1165 end do 1166 1167 ikg = ikg + npw_k 1168 1169 end do ! close loop over ikpt 1170 1171 ipwnsfac = 1 1172 1173 end do ! close loop over ifor 1174 1175 end if ! rfdir(idir) == 1 1176 1177 end do ! close loop over idir 1178 1179 1180 call timab(1008,2,tsec) 1181 call timab(1009,1,tsec) 1182 1183 !Build mpi_enreg%kptdstrb 1184 !array required to communicate the WFs between cpus in berryphase_new.f 1185 !(MPI // over k-points) 1186 if (nproc>1) then 1187 do idir = 1, 3 1188 if (dtset%rfdir(idir) == 1) then 1189 do ifor = 1, 2 1190 1191 ikpt_loc = 0 1192 do isppol = 1, nsppol 1193 1194 do ikpt = 1, dtefield%fnkpt 1195 1196 ikpti = dtefield%indkk_f2ibz(ikpt,1) 1197 nband_k = dtset%nband(ikpti) 1198 ikpt1f = dtefield%ikpt_dk(ikpt,ifor,idir) 1199 ikpt1i = dtefield%indkk_f2ibz(ikpt1f,1) 1200 1201 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpti,1,nband_k,isppol,me)) cycle 1202 1203 ikpt_loc = ikpt_loc + 1 1204 mpi_enreg%kptdstrb(me + 1,ifor+2*(idir-1),ikpt_loc) = & 1205 & ikpt1i + (isppol - 1)*nkpt 1206 1207 mpi_enreg%kptdstrb(me+1,ifor+2*(idir-1),& 1208 & ikpt_loc+dtefield%fmkmem_max*nsppol) = & 1209 & ikpt1f + (isppol - 1)*dtefield%fnkpt 1210 1211 end do ! ikpt 1212 end do ! isppol 1213 end do ! ifor 1214 end if ! dtset%rfdir(idir) == 1 1215 end do ! idir 1216 end if ! nproc>1 1217 1218 !build mpi_enreg%kpt_loc2fbz_sp 1219 ikpt_loc = 0 1220 do isppol = 1, nsppol 1221 do ikpt = 1, dtefield%fnkpt 1222 1223 ikpti = dtefield%indkk_f2ibz(ikpt,1) 1224 nband_k = dtset%nband(ikpti) 1225 1226 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpti,1,nband_k,isppol,me)) cycle 1227 1228 ikpt_loc = ikpt_loc + 1 1229 1230 mpi_enreg%kpt_loc2fbz_sp(me, ikpt_loc, 1) = ikpt 1231 mpi_enreg%kpt_loc2fbz_sp(me, ikpt_loc, 2) = isppol 1232 1233 end do 1234 end do 1235 1236 1237 !parallel case only : 1238 !build mpi_enreg%kpt_loc2ibz_sp, dtefield%cgqindex and dtefield%nneigh 1239 if ((fieldflag).and.(nproc>1)) then 1240 ikpt_loc = 0 1241 do isppol = 1, nsppol 1242 do ikpt = 1, nkpt 1243 1244 ikptf = dtefield%i2fbz(ikpt) 1245 nband_k = dtset%nband(ikpti) 1246 1247 neigh(:) = 0 ; icg = 0 ; ikg = 0 ; flag_kpt = 0; icprj = 0 1248 do idir=1, 3 1249 1250 ! skip idir values for which efield_dot(idir) = 0 1251 if (abs(dtefield%efield_dot(idir)) < tol12) cycle 1252 1253 do ifor = 1, 2 1254 1255 flag = 0 1256 1257 ikpt1f = dtefield%ikpt_dk(ikptf,ifor,idir) 1258 ikpt1i = dtefield%indkk_f2ibz(ikpt1f,1) 1259 1260 dtefield%cgqindex(3,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = ikg 1261 ikg = ikg + npwarr(ikpt1i) 1262 1263 ! check if this neighbour is also a previous neighbour 1264 do ineigh = 1, (ifor+2*(idir-1)) 1265 if (neigh(ineigh) == ikpt1i) then 1266 flag = 1 1267 dtefield%cgqindex(1,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = ineigh 1268 dtefield%cgqindex(2,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = & 1269 & dtefield%cgqindex(2,ineigh,ikpt+(isppol-1)*nkpt) 1270 exit 1271 end if 1272 end do 1273 ! create the cgqindex of the neighbour if necessary 1274 if (flag == 0) then 1275 neigh(ifor+2*(idir-1)) = ikpt1i 1276 dtefield%cgqindex(1,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = & 1277 & ifor+2*(idir-1) 1278 dtefield%cgqindex(2,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = icg 1279 if (isppol == 1) dtefield%nneigh(ikpt) = dtefield%nneigh(ikpt) + 1 1280 icg = icg + npwarr(ikpt1i)*dtefield%nspinor*nband_k 1281 end if 1282 end do !ifor 1283 end do !idir 1284 1285 if (.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me))) then 1286 ! ikpt is one of my kpt_loc 1287 ikpt_loc = ikpt_loc + 1 1288 mpi_enreg%kpt_loc2ibz_sp(me, ikpt_loc, 1) = ikpt 1289 mpi_enreg%kpt_loc2ibz_sp(me, ikpt_loc, 2) = isppol 1290 end if 1291 1292 end do !ikpt 1293 end do !isppol 1294 end if !nproc>1 1295 1296 !should be temporary 1297 !unassigned mpi_enreg%kpt_loc2fbz_sp are empty ; inform other cpu (there are better ways...) 1298 mpi_enreg%mkmem(me) = mkmem 1299 !do ii=ikpt_loc+1,dtefield%fmkmem_max 1300 !mpi_enreg%kpt_loc2fbz_sp(me, ii, 1) = -1 1301 !end do 1302 1303 1304 !(same as mpi_enreg%kptdstrb but for k-points in the iBZ), 1305 !dtefield%cgqindex and dtefield%nneigh 1306 1307 if ((fieldflag).and.(nproc>1)) then 1308 1309 ikpt_loc = 1 1310 do isppol = 1, nsppol 1311 do ikpt = 1, nkpt 1312 1313 nband_k = dtset%nband(ikpt) 1314 ikptf = dtefield%i2fbz(ikpt) 1315 1316 neigh(:) = 0 ; icg = 0 ; ikg = 0 ; flag_kpt = 0; icprj = 0 1317 do idir = 1, 3 1318 1319 ! Skip idir values for which efield_dot(idir) = 0 1320 if (abs(dtefield%efield_dot(idir)) < tol12 .and. (fieldflag)) cycle 1321 1322 do ifor = 1, 2 1323 1324 ikpt1f = dtefield%ikpt_dk(ikptf,ifor,idir) 1325 ikpt1i = dtefield%indkk_f2ibz(ikpt1f,1) 1326 1327 ! dtefield%cgqindex(3,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = ikg 1328 ikg = ikg + npwarr(ikpt1i) 1329 1330 flag = 0 1331 do ineigh = 1, (ifor+2*(idir-1)) 1332 if (neigh(ineigh) == ikpt1i) then 1333 flag = 1 1334 ! dtefield%cgqindex(1,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = ineigh 1335 ! dtefield%cgqindex(2,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = & 1336 ! & dtefield%cgqindex(2,ineigh,ikpt+(isppol-1)*nkpt) 1337 exit 1338 end if 1339 end do 1340 if (flag == 0) then 1341 ! neigh(ifor+2*(idir-1)) = ikpt1i 1342 ! dtefield%cgqindex(1,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = & 1343 ! & ifor+2*(idir-1) 1344 ! dtefield%cgqindex(2,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = icg 1345 ! if (isppol == 1) dtefield%nneigh(ikpt) = dtefield%nneigh(ikpt) + 1 1346 ! icg = icg + npwarr(ikpt1i)*dtset%nspinor*nband_k 1347 end if 1348 1349 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) cycle 1350 1351 flag_kpt = 1 1352 1353 ! MVeithen: the if condition allows to avoid that the same wavefunction 1354 ! is send several times to a particular cpu 1355 1356 end do ! ifor 1357 end do ! idir 1358 1359 if (flag_kpt == 1) ikpt_loc = ikpt_loc + 1 1360 1361 end do ! ikpt 1362 end do ! isppol 1363 1364 end if ! fieldflag 1365 1366 call xmpi_sum(mpi_enreg%kptdstrb,spaceComm,ierr) 1367 call xmpi_sum(mpi_enreg%kpt_loc2fbz_sp,spaceComm,ierr) 1368 if (fieldflag) then 1369 call xmpi_sum(mpi_enreg%kpt_loc2ibz_sp,spaceComm,ierr) 1370 call xmpi_sum(mpi_enreg%mkmem,spaceComm,ierr) 1371 end if 1372 1373 !------------------------------------------------------------------------------ 1374 !------------------------ Estimate critical field ----------------------------- 1375 !------------------------------------------------------------------------------ 1376 1377 !Compute the minimal value of the bandgap required to be below 1378 !the critical field as defined by the relation 1379 !| E_i*a_i | < E_g/n_i 1380 1381 if (fieldflag) then 1382 1383 do idir = 1, 3 1384 ! eg_dir(idir) = abs(dtefield%efield_dot(idir))*dtefield%nkstr(idir) 1385 eg_dir(idir) = abs(dtset%red_efieldbar(idir))*dtefield%nkstr(idir) 1386 end do 1387 1388 1389 eg = maxval(eg_dir) 1390 eg_ev = eg*Ha_eV 1391 1392 if (dtset%optcell ==0 .and. (dtset%berryopt == 4 .or. dtset%berryopt == 14)) then 1393 write(message,'(a,a,a,a,a,a,a,a,f7.2,a,a)')ch10,& 1394 & ' initberry: COMMENT - ',ch10,& 1395 & ' As a rough estimate,',ch10,& 1396 & ' to be below the critical field, the bandgap of your system',ch10,& 1397 & ' should be larger than ',eg_ev,' eV.',ch10 1398 call wrtout(ab_out,message,'COLL') 1399 call wrtout(std_out,message,'COLL') 1400 1401 else 1402 1403 write(message,'(a,a,a,a,a,a,a)') ch10,& 1404 & ' initberry: COMMENT - ',ch10,& 1405 & ' The estimation of critical electric field should be checked after calculation.',ch10,& 1406 & ' It is printed out just after total energy.' ,ch10 1407 1408 call wrtout(ab_out,message,'COLL') 1409 call wrtout(std_out,message,'COLL') 1410 1411 end if 1412 1413 end if 1414 1415 ABI_DEALLOCATE(kg1_k) 1416 1417 call timab(1009,2,tsec) 1418 call timab(1001,2,tsec) 1419 1420 DBG_EXIT("COLL") 1421 1422 end subroutine initberry