TABLE OF CONTENTS
ABINIT/elpolariz [ Functions ]
NAME
elpolariz
FUNCTION
Calculate corrections to total energy from polarising electric field with or without Berry phases (berryopt keyword)
INPUTS
atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f) cg(2,mcg)=planewave coefficients of wavefunctions cprj(natom,mcprj*usecprj)=<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 gprimd(3,3)=reciprocal space dimensional primitive translations hdr <type(hdr_type)>=the header of wf, den and pot files kg(3,mpw*mkmem)=reduced planewave coordinates mband=maximum number of bands mcg=size of wave-functions array (cg) =mpw*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 nattyp(ntypat)= # atoms of each type. 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 option = 1: compute Berryphase polarization 2: compute finite difference expression of the ddk 3: compute polarization & ddk pawrhoij(my_natom*usepaw) <type(pawrhoij_type)> atomic occupancies pawtab(dtset%ntypat) <type(pawtab_type)>=paw tabulated starting data pel_cg(3) = reduced coordinates of the electronic polarization (a. u.) computed in the SCF loop pelev(3)= expectation value polarization term (PAW only) in cartesian coordinates 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 rprimd(3,3)=dimensional real space primitive translations (bohr) ucvol=unit cell volume in bohr**3. usecprj=1 if cprj datastructure has been allocated xred(3,natom)=reduced atomic coordinates
OUTPUT
(see side effects)
SIDE EFFECTS
dtefield <type(efield_type)> = variables related to Berry phase and electric field calculations (see initberry.f). In case berryopt = 4/6/7/14/16/17, the overlap matrices computed in this routine are stored in dtefield%smat in order to be used in the electric field calculation. enefield=field energy etotal=total energy, might be correct by improved polarization computation pel(3) = reduced coordinates of the electronic polarization (a. u.) pion(3)= reduced coordinates of the ionic polarization (a. u.)
SOURCE
122 subroutine elpolariz(atindx1,cg,cprj,dtefield,dtfil,dtset,etotal,enefield,gprimd,hdr,& 123 & kg,mband,mcg,mcprj,mkmem,mpi_enreg,mpw,my_natom,natom,nattyp,nkpt,& 124 & npwarr,nsppol,ntypat,pawrhoij,pawtab,& 125 & pel,pel_cg,pelev,pion,psps,pwind,pwind_alloc,& 126 & pwnsfac,rprimd,ucvol,usecprj,xred) 127 128 !Arguments ------------------------------------ 129 !scalars 130 integer,intent(in) :: mband,mcg,mcprj,mkmem,mpw,my_natom,natom,nkpt,nsppol,ntypat 131 integer,intent(in) :: pwind_alloc,usecprj 132 real(dp),intent(in) :: ucvol 133 real(dp),intent(inout) :: enefield,etotal 134 type(MPI_type),intent(in) :: mpi_enreg 135 type(datafiles_type),intent(in) :: dtfil 136 type(dataset_type),intent(inout) :: dtset 137 type(efield_type),intent(inout) :: dtefield 138 type(hdr_type),intent(inout) :: hdr 139 type(pseudopotential_type),intent(in) :: psps 140 !arrays 141 integer,intent(in) :: atindx1(natom),kg(3,mpw*mkmem),nattyp(ntypat) 142 integer,intent(in) :: npwarr(nkpt),pwind(pwind_alloc,2,3) 143 real(dp),intent(in) :: cg(2,mcg),gprimd(3,3) 144 real(dp),intent(in) :: pel_cg(3),pwnsfac(2,pwind_alloc),rprimd(3,3) 145 real(dp),intent(inout) :: pel(3),pelev(3),pion(3),xred(3,natom) 146 type(pawcprj_type),intent(in) :: cprj(natom,mcprj*usecprj) 147 type(pawrhoij_type), intent(in) :: pawrhoij(my_natom*psps%usepaw) 148 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*psps%usepaw) 149 150 !Local variables------------------------------- 151 !scalars 152 integer :: mcg1_3,my_nspinor,option,unit_out,iir,jjr,kkr 153 real(dp) :: pdif_mod,eenth,ucvol_local 154 logical :: save_cg1_3 155 character(len=500) :: message 156 !arrays 157 real(dp) :: gmet(3,3),gprimdlc(3,3),pdif(3),ptot(3),red_ptot(3),rmet(3,3) 158 !! ptot(3) = total polarization (not reduced) REC 159 !! red_ptot(3) = internal reduced total polarization REC 160 real(dp) :: A(3,3),A1(3,3),A_new(3,3),efield_new(3) 161 real(dp),allocatable :: cg1_3(:,:,:) 162 163 ! ************************************************************************* 164 165 DBG_ENTER("COLL") 166 167 if (usecprj==0.and.psps%usepaw==1) then 168 write (message,'(3a)')& 169 & 'cprj datastructure must be allocated !',ch10,& 170 & 'Action: change pawusecp input keyword.' 171 ABI_ERROR(message) 172 end if 173 174 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 175 176 if(dtset%berryopt>0 .and. dtset%berryopt/=4 .and. dtset%berryopt/=6 .and. dtset%berryopt/=7 .and. & 177 & dtset%berryopt/=14 .and. dtset%berryopt/=16 .and. dtset%berryopt/=17)then !!HONG 178 179 if (dtset%berryopt==1 .or. dtset%berryopt==3) then 180 call berryphase(atindx1,dtset%bdberry,cg,gprimd,dtset%istwfk,& 181 & dtset%kberry,kg,dtset%kptns,dtset%kptopt,dtset%kptrlatt,& 182 & mband,mcg,mkmem,mpw,natom,nattyp,dtset%nband,dtset%nberry,npwarr,& 183 & my_nspinor,nsppol,psps%ntypat,nkpt,rprimd,ucvol,& 184 & xred,psps%ziontypat) 185 end if 186 187 if (dtset%berryopt==2 .or. dtset%berryopt==3) then 188 call uderiv(dtset%bdberry,cg,gprimd,hdr,dtset%istwfk,& 189 & dtset%kberry,kg,dtset%kptns,dtset%kptopt,& 190 & dtset%kptrlatt,mband,mcg,mkmem,mpi_enreg,mpw,& 191 & natom,dtset%nband,dtset%nberry,npwarr,my_nspinor,nsppol,& 192 & nkpt,dtfil%unddk,dtfil%fnameabo_1wf) 193 end if 194 195 else if(dtset%berryopt<0 .or. dtset%berryopt==4 .or. dtset%berryopt==6 .or. dtset%berryopt==7 .or. & 196 & dtset%berryopt==14 .or. dtset%berryopt==16 .or. dtset%berryopt==17)then !!HONG 197 198 select case (dtset%berryopt) 199 case (-5) 200 option = 2 201 case (-3) 202 option = 3 203 case (-2) 204 option = 2 205 case (-1) 206 option = 1 207 case (4) 208 option = 1 209 pel(:) = zero 210 pelev(:) = zero 211 case (6) !!HONG 212 option = 1 213 pel(:) = zero 214 pelev(:) = zero 215 case (7) !!HONG 216 option = 1 217 pel(:) = zero 218 pelev(:) = zero 219 case (14) !!HONG 220 option = 1 221 pel(:) = zero 222 pelev(:) = zero 223 case (16) !!HONG 224 option = 1 225 pel(:) = zero 226 pelev(:) = zero 227 case (17) !!HONG 228 option = 1 229 pel(:) = zero 230 pelev(:) = zero 231 end select 232 233 unit_out = ab_out 234 235 ! by default berryphase_new writes ddk wavefunctions to disk 236 ! with save_cg1_3 set to true it can return them in memory 237 ! this feature is used at present when orbmag is also called, 238 ! as in afterscfloop 239 save_cg1_3 = .FALSE. 240 mcg1_3 = 0 241 ABI_MALLOC(cg1_3,(2,mcg1_3,3)) 242 call berryphase_new(atindx1,cg,cg1_3,cprj,dtefield,dtfil,dtset,psps,& 243 & gprimd,hdr,psps%indlmn,kg,& 244 & psps%lmnmax,mband,mcg,mcg1_3,mcprj,mkmem,mpi_enreg,mpw,my_natom,natom,npwarr,& 245 & nsppol,psps%ntypat,nkpt,option,pawrhoij,& 246 & pawtab,pel,pelev,pion,ptot,red_ptot,pwind,& !!REC 247 & pwind_alloc,pwnsfac,rprimd,save_cg1_3,dtset%typat,ucvol,& 248 & unit_out,usecprj,psps%usepaw,xred,psps%ziontypat) 249 ABI_FREE(cg1_3) 250 251 dtefield%red_ptot1(:)=red_ptot(:) 252 253 if (dtset%berryopt == 4 .or. dtset%berryopt == 6 .or. dtset%berryopt == 7 .or. & 254 & dtset%berryopt == 14 .or. dtset%berryopt == 16 .or. dtset%berryopt == 17 ) then !!HONG 255 256 ! Check if pel has the same value as pel_cg 257 ! if (psps%usepaw == 1) pel(:) = pel(:) + pelev(:) ! add on-site term for PAW 258 ! if (psps%usepaw == 1) red_ptot(:) = red_ptot(:) + pelev(:) ! add on-site term for PAW !! REC 259 ! above line suppressed because in the PAW case, pel already includes all on-site 260 ! terms and pelev should not be added in additionally. We are computing pelev separately for 261 ! reporting purposes only. 262 ! 13 June 2012 J Zwanziger 263 264 pdif(:) = pel_cg(:) - pel(:) 265 pdif_mod = pdif(1)**2 + pdif(2)**2 + pdif(3)**2 266 267 if (pdif_mod > tol8) then 268 write(message,'(11(a),e16.9)')ch10,& 269 & ' scfcv (electric field calculation) : WARNING -',ch10,& 270 & ' The difference between pel (electronic Berry phase updated ',ch10,& 271 & ' at each SCF cycle)',ch10,& 272 & ' and pel_cg (electronic Berryphase computed using the ',& 273 & 'berryphase routine) is',ch10,& 274 & ' pdif_mod = ',pdif_mod 275 call wrtout(std_out,message,'COLL') 276 write(message,'(a,6(a,e16.9,a))') ch10,& 277 & 'pel_cg(1) = ',pel_cg(1),ch10,& 278 & 'pel_cg(2) = ',pel_cg(2),ch10,& 279 & 'pel_cg(3) = ',pel_cg(3),ch10,& 280 & 'pel(1) = ',pel(1),ch10,& 281 & 'pel(2) = ',pel(2),ch10,& 282 & 'pel(3) = ',pel(3),ch10 283 ABI_ERROR(message) 284 end if 285 286 ! Use this (more accurate) value of P to recompute enefield 287 if (dtset%berryopt == 4 .or. dtset%berryopt == 14 ) then !!HONG 288 etotal = etotal - enefield 289 290 enefield = -dot_product(dtset%red_efieldbar,red_ptot) 291 call metric(gmet,gprimdlc,-1,rmet,rprimd,ucvol_local) 292 eenth = zero 293 do iir=1,3 294 do jjr=1,3 295 eenth= eenth+gmet(iir,jjr)*dtset%red_efieldbar(iir)*dtset%red_efieldbar(jjr) !! HONG g^{-1})_ij ebar_i ebar_j 296 end do 297 end do 298 eenth=-1_dp*(ucvol_local/(8.0d0*pi))*eenth 299 enefield=enefield+eenth 300 301 etotal = etotal + enefield 302 303 write(message,'(a,a)')ch10,& 304 & ' Stress tensor under a constant electric field:' 305 call wrtout(std_out,message,'COLL') 306 call wrtout(ab_out,message,'COLL') 307 308 end if 309 310 ! ! In finite D-field case, turn it into internal energy !!HONG 311 if (dtset%berryopt == 6 .or. dtset%berryopt == 16 ) then 312 etotal = etotal - enefield 313 314 enefield=zero 315 call metric(gmet,gprimdlc,-1,rmet,rprimd,ucvol_local) 316 do iir=1,3 317 do jjr=1,3 318 enefield= enefield+gmet(iir,jjr)*dtset%red_efieldbar(iir)*dtset%red_efieldbar(jjr) !! HONG g^{-1})_ij ebar_i ebar_j 319 end do 320 end do 321 enefield= ucvol_local/(8.0d0*pi)*enefield 322 323 etotal = etotal + enefield 324 325 write(message,'(a,a)')ch10,& 326 & ' Stress tensor under a constant electric displacement field:' 327 call wrtout(std_out,message,'COLL') 328 call wrtout(ab_out,message,'COLL') 329 330 end if 331 332 ! HONG calculate internal energy and electric enthalpy for mixed BC case. 333 if ( dtset%berryopt == 17 ) then 334 etotal = etotal - enefield 335 enefield = zero 336 337 call metric(gmet,gprimdlc,-1,rmet,rprimd,ucvol_local) 338 A(:,:)=(4*pi/ucvol_local)*rmet(:,:) 339 A1(:,:)=A(:,:) 340 A_new(:,:)=A(:,:) 341 efield_new(:)=dtset%red_efield(:) 342 eenth = zero 343 344 do kkr=1,3 345 if (dtset%jfielddir(kkr)==1) then ! fixed ebar direction 346 ! step 1 add -ebar*p 347 eenth=eenth - dtset%red_efieldbar(kkr)*red_ptot(kkr) 348 349 ! step 2 chang to e_new (change e to ebar) 350 efield_new(kkr)=dtset%red_efieldbar(kkr) 351 352 ! step 3 chang matrix A to A1 353 354 do iir=1,3 355 do jjr=1,3 356 if (iir==kkr .and. jjr==kkr) A1(iir,jjr)=-1.0/A(kkr,kkr) 357 if ((iir==kkr .and. jjr/=kkr) .or. (iir/=kkr .and. jjr==kkr)) & 358 & A1(iir,jjr)=-1.0*A(iir,jjr)/A(kkr,kkr) 359 if (iir/=kkr .and. jjr/=kkr) A1(iir,jjr)=A(iir,jjr)-A(iir,kkr)*A(kkr,jjr)/A(kkr,kkr) 360 end do 361 end do 362 363 A(:,:)=A1(:,:) 364 A_new(:,:)=A1(:,:) 365 end if 366 367 end do ! end fo kkr 368 369 370 do iir=1,3 371 do jjr=1,3 372 eenth= eenth+(1/2.0)*A_new(iir,jjr)*efield_new(iir)*efield_new(jjr) 373 end do 374 end do 375 376 enefield=eenth 377 etotal = etotal + enefield 378 379 write(message,'(a,a)')ch10,& 380 & ' Stress tensor under a constant (mixed) electric and electric displacement field:' 381 call wrtout(std_out,message,'COLL') 382 call wrtout(ab_out,message,'COLL') 383 384 end if ! berryopt==17 385 386 387 ! MVeithen: to clarify 388 ! Which stress tensor should be used in structural optimizations? 389 ! The one at constant electric field or at constant potential drop. 390 ! write(message,'(a,a)')ch10,& 391 ! & ' Stress tensor imposing a constant electric field:' 392 ! call wrtout(std_out,message,'COLL') 393 ! call wrtout(ab_out,message,'COLL') 394 395 end if ! dtset%berryopt == 4/6/7/14/16/17 396 397 end if ! dtset%berryopt>0 or dtset%berryopt/=4/6/7/14/16/17 398 399 DBG_EXIT("COLL") 400 401 end subroutine elpolariz
ABINIT/m_elpolariz [ Modules ]
NAME
m_elpolariz
FUNCTION
COPYRIGHT
Copyright (C) 2005-2024 ABINIT group (XG, NSAI, MKV) 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
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_elpolariz 23 24 use defs_basis 25 use m_abicore 26 use m_errors 27 use m_efield 28 use m_xmpi 29 use m_wffile 30 use m_hdr 31 use m_dtset 32 use m_dtfil 33 34 use defs_datatypes, only : pseudopotential_type 35 use defs_abitypes, only : MPI_type 36 use m_geometry, only : metric 37 use m_symtk, only : matr3inv 38 use m_hide_lapack, only : dzgedi, dzgefa 39 use m_rwwf, only : rwwf 40 use m_pawtab, only : pawtab_type 41 use m_pawrhoij, only : pawrhoij_type 42 use m_pawcprj, only : pawcprj_type 43 use m_berryphase, only : berryphase 44 use m_berryphase_new, only : berryphase_new 45 46 implicit none 47 48 private
ABINIT/uderiv [ Functions ]
NAME
uderiv
FUNCTION
This routine is called computes the derivative of ground-state wavefunctions with respect to k (du/dk) by finite differencing on neighbouring k points Work for nsppol=1 or 2, but only accept nspinor=1,
INPUTS
bdberry(4)=band limits for Berry phase contributions (or du/dk) spin up and spin down (bdberry(3:4) is irrelevant when nsppol=1) cg(2,mcg)=planewave coefficients of wavefunctions gprimd(3,3)=reciprocal space dimensional primitive translations hdr <type(hdr_type)>=the header of wf, den and pot files istwfk(nkpt_)=input option parameter that describes the storage of wfs kberry(3,20)= different delta k for Berry phases(or du/dk), in unit of kptrlatt only kberry(1:3,1:nberry) is relevant kg(3,mpw*mkmem)=reduced planewave coordinates kpt_(3,nkpt_)=reduced coordinates of k points generated by ABINIT, kpt_ samples half the BZ if time-reversal symetrie is used kptopt=2 when time-reversal symmetry is used kptrlatt(3,3)=k-point lattice specification mband=maximum number of bands mcg=size of wave-functions array (cg) =mpw*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 natom=number of atoms in cell nband(nkpt*nsppol)=number of bands at each k point, for each polarization nberry=number of Berry phases(or du/dk) to be computed nkpt=number of k points npwarr(nkpt)=number of planewaves in basis at this k point nspinor=number of spinorial components of the wavefunctions (on current proc) nsppol=1 for unpolarized, 2 for spin-polarized unddk=unit number for ddk file
OUTPUT
(the ddk wavefunctions are written on disk)
SIDE EFFECTS
TODO
Cleaning, checking for rules Should allow for time-reversal symmetry (istwfk) WARNING : the use of nspinor is completely erroneous
NOTES
Local Variables:
cmatrix(:,:,:)= overlap matrix of size maxband*maxband
cg_index(:,:,:)= unpacked cg index array for specific band,
k point and polarization.
det(2,2)= intermediate output of Lapack routine zgedi.f
dk(3)= step taken to the next k mesh point along the kberry direction
gpard(3)= dimensionalreciprocal lattice vector G along which the
polarization is computed
kg_kpt(:,:,:)= unpacked reduced planewave coordinates with subscript of
planewave and k point
kpt(3,nkpt)=reduced coordinates of k-point grid that samples the whole BZ
kpt_flag(nkpt)=kpt_flag(ikpt)=0 when the wf was generated by the ABINIT
code
kpt_flag(ikpt) gives the indices of the k-point
related to ikpt by time revers
maxband/minband= control the minimum and maximum band calculated in the
overlap matrix
npw_k= npwarr(ikpt), number of planewaves in basis at this k point
shift_g_2(nkpt,nkpt)= .true. if the k point should be shifted by a G vector;
.false. if not
tr(2)=variable that changes k to -k
G to -G
$c_g$ to $c_g^*$ when time-reversal symetrie is used
SOURCE
479 subroutine uderiv(bdberry,cg,gprimd,hdr,istwfk,kberry,kg,kpt_,kptopt,kptrlatt,& 480 & mband,mcg,mkmem,mpi_enreg,mpw,natom,nband,nberry,npwarr,nspinor,nsppol,nkpt_,& 481 & unddk,fnameabo_1wf) 482 483 !Arguments ------------------------------------ 484 !scalars 485 integer,intent(in) :: kptopt,mband,mcg,mkmem,mpw,natom,nberry,nkpt_,nspinor 486 integer,intent(in) :: nsppol,unddk 487 type(MPI_type),intent(in) :: mpi_enreg 488 type(hdr_type),intent(inout) :: hdr 489 !arrays 490 integer,intent(in) :: bdberry(4),istwfk(nkpt_),kberry(3,20),kg(3,mpw*mkmem) 491 integer,intent(in) :: kptrlatt(3,3),nband(nkpt_*nsppol),npwarr(nkpt_) 492 real(dp),intent(in) :: cg(2,mcg),gprimd(1:3,1:3) 493 real(dp),intent(in) :: kpt_(3,nkpt_) 494 character(len=fnlen),intent(in) :: fnameabo_1wf 495 496 !Local variables ------------------------- 497 !scalars 498 integer,parameter :: master=0 499 integer :: iomode,band_in,cg_index_iband,fform,flag1 500 integer :: formeig,iband,iberry,icg,idir,ierr,ifor,ii,ikpt,ikpt2,ikpt_ 501 integer :: index,index1,info,ipert,ipw,isgn,isppol,jband,jj,jkpt,jkpt_ 502 integer :: maxband,mcg_disk,me,minband,nband_diff,nband_k 503 integer :: nkpt,npw_k,pertcase,rdwr,read_k,spaceComm 504 integer :: tim_rwwf 505 real(dp) :: gmod,twodk 506 character(len=500) :: message 507 character(len=fnlen) :: fiwf1o 508 type(wffile_type) :: wffddk 509 !arrays 510 integer :: kg_jl(0,0,0),kpt_flag(2*nkpt_) 511 integer,allocatable :: cg_index(:,:,:),ikpt_dk(:,:),ipvt(:) 512 integer,allocatable :: kg_kpt(:,:,:) 513 real(dp) :: det(2,2),diffk(3),diffk2(3),dk(3),gpard(3),klattice(3,3) 514 real(dp) :: kptrlattr(3,3),tr(2) 515 real(dp) :: cg_disk(0,0,0) 516 real(dp),allocatable :: cmatrix(:,:,:),dudk(:,:) 517 real(dp),allocatable :: eig_dum_2(:),kpt(:,:) 518 real(dp),allocatable :: occ_dum_2(:),phi(:,:,:),u_tilde(:,:,:,:),zgwork(:,:) 519 logical,allocatable :: shift_g_2(:,:) 520 521 ! ************************************************************************* 522 523 if(min(2,(1+mpi_enreg%paral_spinor)*nspinor)==2)then 524 ABI_ERROR('uderiv: does not yet work for nspinor=2') 525 end if 526 527 if(maxval(istwfk(:))/=1)then 528 write(message,'(3a)')& 529 & 'Sorry, this routine does not work yet with istwfk/=1.',ch10,& 530 & 'This should have been tested previously ...' 531 ABI_BUG(message) 532 end if 533 534 if (kptopt==3) then 535 nkpt = nkpt_ 536 ABI_MALLOC(kpt,(3,nkpt)) 537 kpt(:,:)=kpt_(:,:) 538 else if (kptopt==2) then 539 nkpt = nkpt_*2 540 ABI_MALLOC(kpt,(3,nkpt)) 541 do ikpt = 1,nkpt/2 542 kpt_flag(ikpt) = 0 543 kpt(:,ikpt)=kpt_(:,ikpt) 544 end do 545 index = 0 546 do ikpt = (nkpt/2+1),nkpt 547 flag1 = 0 548 do jkpt = 1, nkpt/2 549 if (((abs(kpt_(1,ikpt-nkpt/2)+kpt_(1,jkpt))<1.0d-8).or.& 550 & (abs(1-abs(kpt_(1,ikpt-nkpt/2)+kpt_(1,jkpt)))<1.0d-8))& 551 & .and.((abs(kpt_(2,ikpt-nkpt/2)+kpt_(2,jkpt))<1.0d-8).or.& 552 & (abs(1-abs(kpt_(2,ikpt-nkpt/2)+kpt_(2,jkpt)))<1.0d-8))& 553 & .and.((abs(kpt_(3,ikpt-nkpt/2)+kpt_(3,jkpt))<1.0d-8).or.& 554 & (abs(1-abs(kpt_(3,ikpt-nkpt/2)+kpt_(3,jkpt)))<1.0d-8))) then 555 flag1 = 1 556 index = index + 1 557 exit 558 end if 559 end do 560 if (flag1==0) then 561 kpt_flag(ikpt-index)=ikpt-nkpt/2 562 kpt(:,ikpt-index)=-kpt_(:,ikpt-nkpt/2) 563 end if 564 end do 565 nkpt = nkpt - index 566 end if 567 568 !DEBUG 569 !write(101,*) 'beginning write kpt' 570 !do ikpt=1,nkpt 571 !write(101,*) kpt(:,ikpt) 572 !end do 573 !ENDDEBUG 574 575 ABI_MALLOC(shift_g_2,(nkpt,nkpt)) 576 577 !Compute primitive vectors of the k point lattice 578 !Copy to real(dp) 579 kptrlattr(:,:)=kptrlatt(:,:) 580 !Go to reciprocal space (in reduced coordinates) 581 call matr3inv(kptrlattr,klattice) 582 583 do iberry=1,nberry 584 585 ! ************************************************************************** 586 ! Determine the appended index for ddk 1WF files 587 588 do idir=1,3 589 if (kberry(idir,iberry) ==1) then 590 ipert=natom+1 591 pertcase=idir+(ipert-1)*3 592 end if 593 end do 594 595 ! open ddk 1WF file 596 formeig=1 597 598 call appdig(pertcase,fnameabo_1wf,fiwf1o) 599 !call wfk_open_read(wfk, fiwf1o, formeig, iomode, unddk, spaceComm) 600 601 spaceComm=xmpi_comm_self; me=0 ; iomode=IO_MODE_FORTRAN 602 call WffOpen(iomode,spaceComm,fiwf1o,ierr,wffddk,master,me,unddk) 603 604 rdwr=2 ; fform=2 605 call hdr_io(fform,hdr,rdwr,wffddk) 606 607 ! Define offsets, in case of MPI I/O 608 call xdefineOff(formeig,wffddk,mpi_enreg,nband,npwarr,nspinor,nsppol,nkpt_) 609 610 ! ***************************************************************************** 611 ! Calculate dimensional recip lattice vector along which P is calculated 612 ! dk = step to the nearest k point along that direction 613 ! in reduced coordinates 614 615 dk(:)=dble(kberry(1,iberry))*klattice(:,1)+& 616 & dble(kberry(2,iberry))*klattice(:,2)+& 617 & dble(kberry(3,iberry))*klattice(:,3) 618 619 do idir=1,3 620 if (dk(idir)/=0) then 621 twodk=2*dk(idir) 622 end if 623 end do 624 625 gpard(:)=dk(1)*gprimd(:,1)+dk(2)*gprimd(:,2)+dk(3)*gprimd(:,3) 626 gmod=sqrt(dot_product(gpard,gpard)) 627 628 ! ****************************************************************************** 629 ! Select the k grid points along the direction to compute dudk 630 ! dk = step to the nearest k point along that direction 631 632 ! For each k point, find k_prim such that k_prim= k + dk mod(G) 633 ! where G is a vector of the reciprocal lattice 634 ABI_MALLOC(ikpt_dk,(2,nkpt)) 635 ikpt_dk(1:2,1:nkpt)=0 636 shift_g_2(:,:)= .false. 637 638 do ikpt=1,nkpt 639 do ikpt2=1,nkpt 640 diffk(:)=abs(kpt(:,ikpt2)-kpt(:,ikpt)-dk(:)) 641 diffk2(:)=abs(kpt(:,ikpt2)-kpt(:,ikpt)+dk(:)) 642 if (sum(abs(diffk(:)-nint(diffk(:))))<3*tol8)then 643 ikpt_dk(1,ikpt)=ikpt2 644 if(sum(diffk(:))>=3*tol8) shift_g_2(ikpt,ikpt2) = .true. 645 end if 646 if (sum(abs(diffk2(:)-nint(diffk2(:))))<3*tol8)then 647 ikpt_dk(2,ikpt)=ikpt2 648 if(sum(diffk2(:))>=3*tol8) shift_g_2(ikpt,ikpt2) = .true. 649 end if 650 end do 651 end do 652 653 write(message,'(a,a,a,3f9.5,a,a,3f9.5,a)')ch10,& 654 & ' Computing the derivative for reciprocal vector:',ch10,& 655 & dk(:),' (in reduced coordinates)',ch10,& 656 & gpard(1:3),' (in cartesian coordinates - atomic units)' 657 call wrtout(ab_out,message,'COLL') 658 call wrtout(std_out,message,'COLL') 659 660 if(nsppol==1)then 661 write(message, '(a,i5,a,i5)')& 662 & ' From band number',bdberry(1),' to band number',bdberry(2) 663 else 664 write(message, '(a,i5,a,i5,a,a,a,i5,a,i5,a)')& 665 & ' From band number',bdberry(1),' to band number',bdberry(2),' for spin up,',& 666 & ch10,& 667 & ' from band number',bdberry(3),' to band number',bdberry(4),' for spin down.' 668 end if 669 call wrtout(ab_out,message,'COLL') 670 call wrtout(std_out,message,'COLL') 671 672 ! ***************************************************************************** 673 ABI_MALLOC(dudk,(2,mpw*nspinor*mband*nsppol)) 674 ABI_MALLOC(eig_dum_2,((2*mband)**formeig*mband)) 675 ABI_MALLOC(occ_dum_2,((2*mband)**formeig*mband)) 676 dudk(1:2,:)=0.0_dp 677 eig_dum_2=0.0_dp 678 occ_dum_2=0.0_dp 679 680 if (mkmem/=0) then 681 682 ! Find the location of each wavefunction 683 684 ABI_MALLOC(cg_index,(mband,nkpt_,nsppol)) 685 icg = 0 686 do isppol=1,nsppol 687 do ikpt=1,nkpt_ 688 nband_k=nband(ikpt+(isppol-1)*nkpt_) 689 npw_k=npwarr(ikpt) 690 do iband=1,nband_k 691 cg_index(iband,ikpt,isppol)=(iband-1)*npw_k*nspinor+icg 692 end do 693 icg=icg+npw_k*nspinor*nband(ikpt) 694 end do 695 end do 696 697 ! Find the planewave vectors for each k point 698 ! SHOULD BE REMOVED WHEN ANOTHER INDEXING TECHNIQUE WILL BE USED FOR kg 699 ABI_MALLOC(kg_kpt,(3,mpw*nspinor,nkpt_)) 700 kg_kpt(:,:,:) = 0 701 index1 = 0 702 do ikpt=1,nkpt_ 703 npw_k=npwarr(ikpt) 704 do ipw=1,npw_k*nspinor 705 kg_kpt(1:3,ipw,ikpt)=kg(1:3,ipw+index1) 706 end do 707 index1=index1+npw_k*nspinor 708 end do 709 end if 710 711 ! ************************************************************************* 712 ! Loop over spins 713 do isppol=1,nsppol 714 715 minband=bdberry(2*isppol-1) 716 maxband=bdberry(2*isppol) 717 718 if(minband<1)then 719 write(message,'(a,i0,a)')' The band limit minband= ',minband,', is lower than 0.' 720 ABI_BUG(message) 721 end if 722 723 if(maxband<1)then 724 write(message,'(a,i0,a)')' The band limit maxband= ',maxband,', is lower than 0.' 725 ABI_BUG(message) 726 end if 727 728 if(maxband<minband)then 729 write(message,'(a,i0,a,i0)')' maxband= ',maxband,', is lower than minband= ',minband 730 ABI_BUG(message) 731 end if 732 733 ! Loop over k points 734 do ikpt_=1,nkpt_ 735 736 read_k = 0 737 738 ikpt=ikpt_ 739 tr(1) = 1.0_dp 740 741 if (kptopt==2) then 742 if (read_k == 0) then 743 if (kpt_flag(ikpt_)/=0) then 744 tr(1) = -1.0_dp 745 ikpt= kpt_flag(ikpt_) 746 end if 747 else !read_k 748 if (kpt_flag(ikpt_)/=0) then 749 tr(-1*read_k+3) = -1.0_dp 750 ikpt= kpt_flag(ikpt_) 751 end if 752 end if !read_k 753 end if !kptopt 754 755 nband_k=nband(ikpt+(isppol-1)*nkpt_) 756 757 if(nband_k<maxband)then 758 write(message,'(a,i0,a,i0)')' maxband=',maxband,', is larger than nband(i,isppol)=',nband_k 759 ABI_BUG(message) 760 end if 761 762 npw_k=npwarr(ikpt) 763 764 ABI_MALLOC(u_tilde,(2,npw_k,maxband,2)) 765 u_tilde(1:2,1:npw_k,1:maxband,1:2)=0.0_dp 766 767 ! ifor = 1,2 represents forward and backward neighbouring k points of ikpt 768 ! respectively along dk direction 769 770 do ifor=1,2 771 772 ABI_MALLOC(phi,(2,mpw,mband)) 773 ABI_MALLOC(cmatrix,(2,maxband,maxband)) 774 phi(1:2,1:mpw,1:mband)=0.0_dp; cmatrix(1:2,1:maxband,1:maxband)=0.0_dp 775 776 isgn=(-1)**ifor 777 jkpt_= ikpt_dk(ifor,ikpt_) 778 779 tr(2) = 1.0_dp 780 781 jkpt=jkpt_ 782 783 if (kptopt==2) then 784 if (read_k == 0) then 785 if (kpt_flag(jkpt_)/=0) then 786 tr(2) = -1.0_dp 787 jkpt= kpt_flag(jkpt_) 788 end if 789 else !read_k 790 if (kpt_flag(jkpt_)/=0) then 791 tr(read_k) = -1.0_dp 792 jkpt= kpt_flag(jkpt_) 793 end if 794 end if !read_k 795 end if !kptopt 796 797 if (ifor==1) read_k = 2 798 799 jj = read_k 800 ii = -1*read_k+3 801 802 call waveformat(cg,cg_disk,cg_index,phi,dk,ii,ikpt,& 803 & ikpt_,isgn,isppol,jj,jkpt,jkpt_,kg_kpt,kpt,kg_jl,maxband,mband,mcg,mcg_disk,& 804 & minband,mkmem,mpw,nkpt,nkpt_,npwarr,nsppol,nspinor,shift_g_2,tr) 805 806 ! Compute the overlap matrix <u_k|u_k+b> 807 808 do iband=minband,maxband 809 cg_index_iband=cg_index(iband,ikpt,isppol) 810 do jband=minband,maxband 811 do ipw=1,npwarr(ikpt) 812 cmatrix(1,iband,jband)=cmatrix(1,iband,jband)+& 813 & cg(1,ipw+cg_index_iband)*phi(1,ipw,jband)+& 814 & tr(ii)*cg(2,ipw+cg_index_iband)*tr(jj)*phi(2,ipw,jband) 815 816 cmatrix(2,iband,jband)=cmatrix(2,iband,jband)+& 817 & cg(1,ipw+cg_index_iband)*tr(jj)*phi(2,ipw,jband)-& 818 & tr(ii)*cg(2,ipw+cg_index_iband)*phi(1,ipw,jband) 819 end do 820 end do 821 end do 822 823 ! Compute the inverse of cmatrix(1:2,minband:maxband, minband:maxband) 824 825 band_in = maxband - minband + 1 826 ABI_MALLOC(ipvt,(maxband)) 827 ABI_MALLOC(zgwork,(2,1:maxband)) 828 829 ! Last argument of zgedi means calculate inverse only 830 call dzgefa(cmatrix(1,minband,minband),maxband, band_in,ipvt,info) 831 call dzgedi(cmatrix(1,minband,minband),maxband, band_in,ipvt,det,zgwork,01) 832 833 ABI_FREE(zgwork) 834 ABI_FREE(ipvt) 835 836 ! Compute the product of Inverse overlap matrix with the wavefunction 837 838 do iband=minband,maxband 839 do ipw=1,npwarr(ikpt) 840 u_tilde(1,ipw,iband,ifor)= & 841 & dot_product(cmatrix(1,minband:maxband,iband),& 842 & phi(1,ipw,minband:maxband))-& 843 & dot_product(cmatrix(2,minband:maxband,iband),& 844 & tr(jj)*phi(2,ipw,minband:maxband)) 845 u_tilde(2,ipw,iband,ifor)= & 846 & dot_product(cmatrix(1,minband:maxband,iband),& 847 & tr(jj)*phi(2,ipw,minband:maxband))+& 848 & dot_product(cmatrix(2,minband:maxband,iband),& 849 & phi(1,ipw,minband:maxband)) 850 end do 851 end do 852 ABI_FREE(cmatrix) 853 ABI_FREE(phi) 854 855 end do !ifor 856 857 ! Compute dudk for ikpt 858 859 npw_k=npwarr(ikpt) 860 861 do iband=minband,maxband 862 863 icg=(iband-minband)*npw_k 864 865 dudk(1,1+icg:npw_k+icg)=(u_tilde(1,1:npw_k,iband,1)-& 866 & u_tilde(1,1:npw_k,iband,2))/twodk 867 868 dudk(2,1+icg:npw_k+icg)=(u_tilde(2,1:npw_k,iband,1)-& 869 & u_tilde(2,1:npw_k,iband,2))/twodk 870 871 end do 872 873 tim_rwwf=0 874 mcg_disk=mpw*nspinor*mband 875 nband_diff=maxband-minband+1 876 call rwwf(dudk,eig_dum_2,formeig,0,0,ikpt,isppol,kg_kpt(:,:,ikpt),& 877 & mband,mcg_disk,mpi_enreg,nband_diff,nband_diff,& 878 & npw_k,nspinor,occ_dum_2,2,1,tim_rwwf,wffddk) 879 880 !call wfk_read_band_block(wfk, band_block, ikpt, isppol, sc_mode, 881 ! kg_k=kg_kpt(:,:,ikpt), cg_k=dudk, eig_k=eig_dum, occ_k=occ_dum) 882 883 ABI_FREE(u_tilde) 884 885 end do !ikpt 886 end do !isppol 887 888 ABI_FREE(eig_dum_2) 889 ABI_FREE(occ_dum_2) 890 ABI_FREE(dudk) 891 892 call WffClose(wffddk,ierr) 893 !call wfk_close(wfk) 894 895 ABI_FREE(kg_kpt) 896 ABI_FREE(cg_index) 897 ABI_FREE(ikpt_dk) 898 899 end do ! iberry 900 901 ABI_FREE(shift_g_2) 902 ABI_FREE(kpt) 903 904 write(std_out,*) 'uderiv: exit ' 905 906 end subroutine uderiv
ABINIT/waveformat [ Functions ]
NAME
waveformat
FUNCTION
This routine is to find the matched pairs of plane waves between two neighbouring k points and load a new pw coefficients array cg_new Was written first by Na Sai (thanks), but unfortunately without any comment ...
INPUTS
cg(2,mpw*nspinor*mband*mkmem*nsppol)= input planewave coefficients, in case mkmem/=0
cg_disk(2,mpw*nspinor*mband,2)= input planewave coefficients, in case mkmem==0
cg_index(mband,nkpt_,nsppol)=index of wavefunction iband,ikpt,isppol in the array cg.
dk(3)= step taken to the next k mesh point along the kberry direction (see also isgn)
ii=(to be documented)
ikpt=index of the first k-point in the reduced Brillouin zone
ikpt_=index of the first k-point in the full Brillouin zone
isgn=1 if dk(3) is connecting the k-points (ikpt_ and jkpt)
=-1 if -dk(3) is connecting the k-points
isppol=1 if spin-up, =2 if spin-down
jj=(to be documented)
jkpt=index of the second k-point in the reduced Brillouin zone
jkpt_=index of the second k-point in the full Brillouin zone
kg_kpt(:,:,:)= unpacked reduced planewave coordinates with subscript of
planewave and k point
kpt(3,nkpt)=reduced coordinates of k-point grid that samples the whole BZ
kg_jl(3,mpw,2)=(to be documented)
maxband/minband= control the minimum and maximum band calculated in the
overlap matrix
mband=maximum number of bands (dimension of several cg* arrays)
mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
mcg_disk=size of wave-functions array (cg_disk) =mpw*nspinor*mband
mkmem= if 0, the wavefunctions are input in cg_disk, otherwise in cg
mpw=maximum number of planewaves (dimension of several cg* arrays)
nkpt=number of k points (full Brillouin zone !?!)
nkpt_=number of k points (reduced Brillouin zone !?!)
npwarr(nkpt)=number of planewaves in basis at this k point
nspinor=number of spinorial components of the wavefunctions
nsppol=1 for unpolarized, 2 for spin-polarized
shift_g_2(nkpt,nkpt)=non-zero if a G vector along kberry is needed to connect k points
tr(2)=variable that changes k to -k
G to -G
$c_g$ to $c_g^*$ when time-reversal symetrie is used
OUTPUT
cg_new(2,mpw,maxband)=planewave coefficients transferred onto the set of planewaves at k
SIDE EFFECTS
SOURCE
963 subroutine waveformat(cg,cg_disk,cg_index,cg_new,dk,ii,ikpt,& 964 & ikpt_,isgn,isppol,jj,jkpt,jkpt_,kg_kpt,kpt,kg_jl,maxband,mband,mcg,mcg_disk,& 965 & minband,mkmem,mpw,nkpt,nkpt_,npwarr,nsppol,nspinor,shift_g_2,tr) 966 967 !Arguments ------------------------------------ 968 !scalars 969 integer,intent(in) :: ii,ikpt,ikpt_,isgn,isppol,jj,jkpt,jkpt_,maxband,mband,mcg,mcg_disk 970 integer,intent(in) :: minband,mkmem,mpw,nkpt,nkpt_,nspinor,nsppol 971 !arrays 972 integer,intent(in) :: cg_index(mband,nkpt_,nsppol),kg_jl(3,mpw,2) 973 integer,intent(in) :: kg_kpt(3,mpw*nspinor,nkpt_),npwarr(nkpt_) 974 real(dp),intent(in) :: cg(2,mcg) 975 real(dp),intent(in) :: cg_disk(2,mcg_disk,2),dk(3),kpt(3,nkpt),tr(2) 976 real(dp),intent(out) :: cg_new(2,mpw,maxband) 977 logical,intent(in) :: shift_g_2(nkpt,nkpt) 978 979 !Local variables ------------------------- 980 !scalars 981 integer :: cg_index_iband,iband,ipw,jpw,nomatch,npw_k 982 logical :: found_match 983 !arrays 984 integer :: dg(3) 985 986 ! *********************************************************************** 987 988 npw_k=npwarr(ikpt) 989 990 991 nomatch=0 992 993 !If there is no shift of G-vector between ikpt_ and jkpt_ 994 if(shift_g_2(ikpt_,jkpt_) .eqv. .false.) then 995 996 ! DEBUG 997 ! write(111,*)'pair', ikpt_,jkpt_,'noshift' 998 ! ENDDEBUG 999 1000 ! If the original wavefunction is contained in cg_disk 1001 if(mkmem==0) then 1002 1003 do ipw=1,npw_k 1004 1005 found_match = .false. 1006 1007 do jpw=1,npwarr(jkpt) 1008 if (sum(abs(tr(ii)*kg_jl(:,ipw,ii)-tr(jj)*kg_jl(:,jpw,jj)))<3*tol8)then 1009 do iband=minband, maxband 1010 cg_index_iband=(iband-1)*npwarr(jkpt) 1011 cg_new(1:2,ipw,iband)=cg_disk(1:2,jpw+cg_index_iband,jj) 1012 end do 1013 found_match = .true. 1014 exit 1015 end if 1016 end do 1017 1018 if (found_match .eqv. .false.) then 1019 do iband=minband,maxband 1020 cg_new(1:2,ipw,iband)=zero 1021 end do 1022 nomatch = nomatch + 1 1023 end if 1024 1025 end do 1026 1027 ! Here, the wavefunctions are contained in cg 1028 else 1029 1030 do ipw=1,npw_k 1031 1032 found_match = .false. 1033 1034 do jpw=1,npwarr(jkpt) 1035 if (sum(abs(tr(ii)*kg_kpt(:,ipw,ikpt)-tr(jj)*kg_kpt(:,jpw,jkpt)))<3*tol8)then 1036 do iband=minband, maxband 1037 cg_index_iband=cg_index(iband,jkpt,isppol) 1038 cg_new(1:2,ipw,iband)=cg(1:2,jpw+cg_index_iband) 1039 end do 1040 found_match = .true. 1041 exit 1042 end if 1043 end do 1044 1045 if (found_match .eqv. .false.) then 1046 do iband=minband,maxband 1047 cg_new(1:2,ipw,iband)=(0.0_dp,0.0_dp) 1048 end do 1049 nomatch = nomatch + 1 1050 end if 1051 end do 1052 1053 end if 1054 1055 ! DEBUG 1056 ! write(111,*) 'normal pair nomatch=',nomatch 1057 ! ENDDEBUG 1058 1059 ! If there is a G-vector shift between ikpt_ and jkpt_ 1060 else 1061 1062 ! DEBUG 1063 ! write(111,*) 'pair',ikpt_,jkpt_,' need shift' 1064 ! ENDDEBUG 1065 1066 dg(:) = -1*nint(tr(jj)*kpt(:,jkpt)-tr(ii)*kpt(:,ikpt)+isgn*dk(:)) 1067 1068 ! If the original wavefunction is contained in cg_disk 1069 if(mkmem==0) then 1070 1071 do ipw=1,npw_k 1072 1073 found_match = .false. 1074 1075 do jpw=1,npwarr(jkpt) 1076 if (sum(abs(tr(ii)*kg_jl(:,ipw,ii)-(tr(jj)*kg_jl(:,jpw,jj)-& 1077 & dg(:))))<3*tol8)then 1078 1079 do iband=minband, maxband 1080 cg_index_iband=(iband-1)*npwarr(jkpt) 1081 cg_new(1:2,ipw,iband)=cg_disk(1:2,jpw+cg_index_iband,jj) 1082 end do 1083 found_match = .true. 1084 exit 1085 end if 1086 end do 1087 1088 if (found_match .eqv. .false.) then 1089 do iband=minband,maxband 1090 cg_new(1:2,ipw,iband)=(0.0_dp,0.0_dp) 1091 end do 1092 nomatch = nomatch + 1 1093 end if 1094 end do 1095 1096 ! Here, the wavefunctions are contained in cg 1097 else 1098 1099 do ipw=1,npw_k 1100 1101 found_match = .false. 1102 1103 do jpw=1,npwarr(jkpt) 1104 if (sum(abs(tr(ii)*kg_kpt(:,ipw,ikpt)-(tr(jj)*kg_kpt(:,jpw,jkpt)-& 1105 & dg(:))))<3*tol8)then 1106 do iband=minband, maxband 1107 cg_index_iband=cg_index(iband,jkpt,isppol) 1108 cg_new(1:2,ipw,iband)=cg(1:2,jpw+cg_index_iband) 1109 end do 1110 found_match = .true. 1111 exit 1112 end if 1113 end do 1114 1115 if (found_match .eqv. .false.) then 1116 do iband=minband,maxband 1117 cg_new(1:2,ipw,iband)=zero 1118 end do 1119 nomatch = nomatch + 1 1120 end if 1121 end do 1122 1123 end if 1124 1125 ! DEBUG 1126 ! write(111,*) 'special pair nomatch=',nomatch 1127 ! ENDDEBUG 1128 1129 end if 1130 1131 end subroutine waveformat