TABLE OF CONTENTS
ABINIT/uderiv [ Functions ]
NAME
uderiv
FUNCTION
This routine is called in scfcv.f to compute 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,
COPYRIGHT
Copyright (C) 2001-2018 ABINIT group (NSAI).
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
PARENTS
elpolariz
CHILDREN
appdig,dzgedi,dzgefa,hdr_io,matr3inv,rwwf,waveformat,wffclose,wffopen wrtout,xdefineoff
SOURCE
87 #if defined HAVE_CONFIG_H 88 #include "config.h" 89 #endif 90 91 #include "abi_common.h" 92 93 subroutine uderiv(bdberry,cg,gprimd,hdr,istwfk,kberry,kg,kpt_,kptopt,kptrlatt,& 94 & mband,mcg,mkmem,mpi_enreg,mpw,natom,nband,nberry,npwarr,nspinor,nsppol,nkpt_,& 95 & unddk,fnameabo_1wf) 96 97 use defs_basis 98 use defs_abitypes 99 use m_errors 100 use m_xmpi 101 use m_wffile 102 use m_errors 103 use m_profiling_abi 104 use m_hdr 105 106 use m_abilasi, only : dzgedi, dzgefa 107 108 !This section has been created automatically by the script Abilint (TD). 109 !Do not modify the following lines by hand. 110 #undef ABI_FUNC 111 #define ABI_FUNC 'uderiv' 112 use interfaces_14_hidewrite 113 use interfaces_32_util 114 use interfaces_56_io_mpi 115 use interfaces_67_common, except_this_one => uderiv 116 !End of the abilint section 117 118 implicit none 119 120 !Arguments ------------------------------------ 121 !scalars 122 integer,intent(in) :: kptopt,mband,mcg,mkmem,mpw,natom,nberry,nkpt_,nspinor 123 integer,intent(in) :: nsppol,unddk 124 type(MPI_type),intent(in) :: mpi_enreg 125 type(hdr_type),intent(inout) :: hdr 126 !arrays 127 integer,intent(in) :: bdberry(4),istwfk(nkpt_),kberry(3,20),kg(3,mpw*mkmem) 128 integer,intent(in) :: kptrlatt(3,3),nband(nkpt_*nsppol),npwarr(nkpt_) 129 real(dp),intent(in) :: cg(2,mcg),gprimd(1:3,1:3) 130 real(dp),intent(in) :: kpt_(3,nkpt_) 131 character(len=fnlen),intent(in) :: fnameabo_1wf 132 133 !Local variables ------------------------- 134 !scalars 135 integer,parameter :: master=0 136 integer :: iomode,band_in,cg_index_iband,fform,flag1 137 integer :: formeig,iband,iberry,icg,idir,ierr,ifor,ii,ikpt,ikpt2,ikpt_ 138 integer :: index,index1,info,ipert,ipw,isgn,isppol,jband,jj,jkpt,jkpt_ 139 integer :: maxband,mcg_disk,me,minband,nband_diff,nband_k 140 integer :: nkpt,npw_k,pertcase,rdwr,read_k,spaceComm 141 integer :: tim_rwwf 142 real(dp) :: gmod,twodk 143 character(len=500) :: message 144 character(len=fnlen) :: fiwf1o 145 type(wffile_type) :: wffddk 146 !arrays 147 integer :: kg_jl(0,0,0),kpt_flag(2*nkpt_) 148 integer,allocatable :: cg_index(:,:,:),ikpt_dk(:,:),ipvt(:) 149 integer,allocatable :: kg_kpt(:,:,:) 150 real(dp) :: det(2,2),diffk(3),diffk2(3),dk(3),gpard(3),klattice(3,3) 151 real(dp) :: kptrlattr(3,3),tr(2) 152 real(dp) :: cg_disk(0,0,0) 153 real(dp),allocatable :: cmatrix(:,:,:),dudk(:,:) 154 real(dp),allocatable :: eig_dum_2(:),kpt(:,:) 155 real(dp),allocatable :: occ_dum_2(:),phi(:,:,:),u_tilde(:,:,:,:),zgwork(:,:) 156 logical,allocatable :: shift_g_2(:,:) 157 158 ! ************************************************************************* 159 160 if(min(2,(1+mpi_enreg%paral_spinor)*nspinor)==2)then 161 MSG_ERROR('uderiv: does not yet work for nspinor=2') 162 end if 163 164 if(maxval(istwfk(:))/=1)then 165 write(message,'(3a)')& 166 & 'Sorry, this routine does not work yet with istwfk/=1.',ch10,& 167 & 'This should have been tested previously ...' 168 MSG_BUG(message) 169 end if 170 171 if (kptopt==3) then 172 nkpt = nkpt_ 173 ABI_ALLOCATE(kpt,(3,nkpt)) 174 kpt(:,:)=kpt_(:,:) 175 else if (kptopt==2) then 176 nkpt = nkpt_*2 177 ABI_ALLOCATE(kpt,(3,nkpt)) 178 do ikpt = 1,nkpt/2 179 kpt_flag(ikpt) = 0 180 kpt(:,ikpt)=kpt_(:,ikpt) 181 end do 182 index = 0 183 do ikpt = (nkpt/2+1),nkpt 184 flag1 = 0 185 do jkpt = 1, nkpt/2 186 if (((abs(kpt_(1,ikpt-nkpt/2)+kpt_(1,jkpt))<1.0d-8).or.& 187 & (abs(1-abs(kpt_(1,ikpt-nkpt/2)+kpt_(1,jkpt)))<1.0d-8))& 188 & .and.((abs(kpt_(2,ikpt-nkpt/2)+kpt_(2,jkpt))<1.0d-8).or.& 189 & (abs(1-abs(kpt_(2,ikpt-nkpt/2)+kpt_(2,jkpt)))<1.0d-8))& 190 & .and.((abs(kpt_(3,ikpt-nkpt/2)+kpt_(3,jkpt))<1.0d-8).or.& 191 & (abs(1-abs(kpt_(3,ikpt-nkpt/2)+kpt_(3,jkpt)))<1.0d-8))) then 192 flag1 = 1 193 index = index + 1 194 exit 195 end if 196 end do 197 if (flag1==0) then 198 kpt_flag(ikpt-index)=ikpt-nkpt/2 199 kpt(:,ikpt-index)=-kpt_(:,ikpt-nkpt/2) 200 end if 201 end do 202 nkpt = nkpt - index 203 end if 204 205 !DEBUG 206 !write(101,*) 'beginning write kpt' 207 !do ikpt=1,nkpt 208 !write(101,*) kpt(:,ikpt) 209 !end do 210 !ENDDEBUG 211 212 ABI_ALLOCATE(shift_g_2,(nkpt,nkpt)) 213 214 !Compute primitive vectors of the k point lattice 215 !Copy to real(dp) 216 kptrlattr(:,:)=kptrlatt(:,:) 217 !Go to reciprocal space (in reduced coordinates) 218 call matr3inv(kptrlattr,klattice) 219 220 do iberry=1,nberry 221 222 ! ************************************************************************** 223 ! Determine the appended index for ddk 1WF files 224 225 do idir=1,3 226 if (kberry(idir,iberry) ==1) then 227 ipert=natom+1 228 pertcase=idir+(ipert-1)*3 229 end if 230 end do 231 232 ! open ddk 1WF file 233 formeig=1 234 235 call appdig(pertcase,fnameabo_1wf,fiwf1o) 236 !call wfk_open_read(wfk, fiwf1o, formeig, iomode, unddk, spaceComm) 237 238 spaceComm=xmpi_comm_self; me=0 ; iomode=IO_MODE_FORTRAN 239 call WffOpen(iomode,spaceComm,fiwf1o,ierr,wffddk,master,me,unddk) 240 241 rdwr=2 ; fform=2 242 call hdr_io(fform,hdr,rdwr,wffddk) 243 244 ! Define offsets, in case of MPI I/O 245 call xdefineOff(formeig,wffddk,mpi_enreg,nband,npwarr,nspinor,nsppol,nkpt_) 246 247 ! ***************************************************************************** 248 ! Calculate dimensional recip lattice vector along which P is calculated 249 ! dk = step to the nearest k point along that direction 250 ! in reduced coordinates 251 252 dk(:)=dble(kberry(1,iberry))*klattice(:,1)+& 253 & dble(kberry(2,iberry))*klattice(:,2)+& 254 & dble(kberry(3,iberry))*klattice(:,3) 255 256 do idir=1,3 257 if (dk(idir)/=0) then 258 twodk=2*dk(idir) 259 end if 260 end do 261 262 gpard(:)=dk(1)*gprimd(:,1)+dk(2)*gprimd(:,2)+dk(3)*gprimd(:,3) 263 gmod=sqrt(dot_product(gpard,gpard)) 264 265 ! ****************************************************************************** 266 ! Select the k grid points along the direction to compute dudk 267 ! dk = step to the nearest k point along that direction 268 269 ! For each k point, find k_prim such that k_prim= k + dk mod(G) 270 ! where G is a vector of the reciprocal lattice 271 ABI_ALLOCATE(ikpt_dk,(2,nkpt)) 272 ikpt_dk(1:2,1:nkpt)=0 273 shift_g_2(:,:)= .false. 274 275 do ikpt=1,nkpt 276 do ikpt2=1,nkpt 277 diffk(:)=abs(kpt(:,ikpt2)-kpt(:,ikpt)-dk(:)) 278 diffk2(:)=abs(kpt(:,ikpt2)-kpt(:,ikpt)+dk(:)) 279 if (sum(abs(diffk(:)-nint(diffk(:))))<3*tol8)then 280 ikpt_dk(1,ikpt)=ikpt2 281 if(sum(diffk(:))>=3*tol8) shift_g_2(ikpt,ikpt2) = .true. 282 end if 283 if (sum(abs(diffk2(:)-nint(diffk2(:))))<3*tol8)then 284 ikpt_dk(2,ikpt)=ikpt2 285 if(sum(diffk2(:))>=3*tol8) shift_g_2(ikpt,ikpt2) = .true. 286 end if 287 end do 288 end do 289 290 write(message,'(a,a,a,3f9.5,a,a,3f9.5,a)')ch10,& 291 & ' Computing the derivative for reciprocal vector:',ch10,& 292 & dk(:),' (in reduced coordinates)',ch10,& 293 & gpard(1:3),' (in cartesian coordinates - atomic units)' 294 call wrtout(ab_out,message,'COLL') 295 call wrtout(std_out,message,'COLL') 296 297 if(nsppol==1)then 298 write(message, '(a,i5,a,i5)')& 299 & ' From band number',bdberry(1),' to band number',bdberry(2) 300 else 301 write(message, '(a,i5,a,i5,a,a,a,i5,a,i5,a)')& 302 & ' From band number',bdberry(1),' to band number',bdberry(2),' for spin up,',& 303 & ch10,& 304 & ' from band number',bdberry(3),' to band number',bdberry(4),' for spin down.' 305 end if 306 call wrtout(ab_out,message,'COLL') 307 call wrtout(std_out,message,'COLL') 308 309 ! ***************************************************************************** 310 ABI_ALLOCATE(dudk,(2,mpw*nspinor*mband*nsppol)) 311 ABI_ALLOCATE(eig_dum_2,((2*mband)**formeig*mband)) 312 ABI_ALLOCATE(occ_dum_2,((2*mband)**formeig*mband)) 313 dudk(1:2,:)=0.0_dp 314 eig_dum_2=0.0_dp 315 occ_dum_2=0.0_dp 316 317 if (mkmem/=0) then 318 319 ! Find the location of each wavefunction 320 321 ABI_ALLOCATE(cg_index,(mband,nkpt_,nsppol)) 322 icg = 0 323 do isppol=1,nsppol 324 do ikpt=1,nkpt_ 325 nband_k=nband(ikpt+(isppol-1)*nkpt_) 326 npw_k=npwarr(ikpt) 327 do iband=1,nband_k 328 cg_index(iband,ikpt,isppol)=(iband-1)*npw_k*nspinor+icg 329 end do 330 icg=icg+npw_k*nspinor*nband(ikpt) 331 end do 332 end do 333 334 ! Find the planewave vectors for each k point 335 ! SHOULD BE REMOVED WHEN ANOTHER INDEXING TECHNIQUE WILL BE USED FOR kg 336 ABI_ALLOCATE(kg_kpt,(3,mpw*nspinor,nkpt_)) 337 kg_kpt(:,:,:) = 0 338 index1 = 0 339 do ikpt=1,nkpt_ 340 npw_k=npwarr(ikpt) 341 do ipw=1,npw_k*nspinor 342 kg_kpt(1:3,ipw,ikpt)=kg(1:3,ipw+index1) 343 end do 344 index1=index1+npw_k*nspinor 345 end do 346 end if 347 348 ! ************************************************************************* 349 ! Loop over spins 350 do isppol=1,nsppol 351 352 minband=bdberry(2*isppol-1) 353 maxband=bdberry(2*isppol) 354 355 if(minband<1)then 356 write(message,'(a,i0,a)')' The band limit minband= ',minband,', is lower than 0.' 357 MSG_BUG(message) 358 end if 359 360 if(maxband<1)then 361 write(message,'(a,i0,a)')' The band limit maxband= ',maxband,', is lower than 0.' 362 MSG_BUG(message) 363 end if 364 365 if(maxband<minband)then 366 write(message,'(a,i0,a,i0)')' maxband= ',maxband,', is lower than minband= ',minband 367 MSG_BUG(message) 368 end if 369 370 ! Loop over k points 371 do ikpt_=1,nkpt_ 372 373 read_k = 0 374 375 ikpt=ikpt_ 376 tr(1) = 1.0_dp 377 378 if (kptopt==2) then 379 if (read_k == 0) then 380 if (kpt_flag(ikpt_)/=0) then 381 tr(1) = -1.0_dp 382 ikpt= kpt_flag(ikpt_) 383 end if 384 else !read_k 385 if (kpt_flag(ikpt_)/=0) then 386 tr(-1*read_k+3) = -1.0_dp 387 ikpt= kpt_flag(ikpt_) 388 end if 389 end if !read_k 390 end if !kptopt 391 392 nband_k=nband(ikpt+(isppol-1)*nkpt_) 393 394 if(nband_k<maxband)then 395 write(message,'(a,i0,a,i0)')' maxband=',maxband,', is larger than nband(i,isppol)=',nband_k 396 MSG_BUG(message) 397 end if 398 399 npw_k=npwarr(ikpt) 400 401 ABI_ALLOCATE(u_tilde,(2,npw_k,maxband,2)) 402 u_tilde(1:2,1:npw_k,1:maxband,1:2)=0.0_dp 403 404 ! ifor = 1,2 represents forward and backward neighbouring k points of ikpt 405 ! respectively along dk direction 406 407 do ifor=1,2 408 409 ABI_ALLOCATE(phi,(2,mpw,mband)) 410 ABI_ALLOCATE(cmatrix,(2,maxband,maxband)) 411 phi(1:2,1:mpw,1:mband)=0.0_dp; cmatrix(1:2,1:maxband,1:maxband)=0.0_dp 412 413 isgn=(-1)**ifor 414 jkpt_= ikpt_dk(ifor,ikpt_) 415 416 tr(2) = 1.0_dp 417 418 jkpt=jkpt_ 419 420 if (kptopt==2) then 421 if (read_k == 0) then 422 if (kpt_flag(jkpt_)/=0) then 423 tr(2) = -1.0_dp 424 jkpt= kpt_flag(jkpt_) 425 end if 426 else !read_k 427 if (kpt_flag(jkpt_)/=0) then 428 tr(read_k) = -1.0_dp 429 jkpt= kpt_flag(jkpt_) 430 end if 431 end if !read_k 432 end if !kptopt 433 434 if (ifor==1) read_k = 2 435 436 jj = read_k 437 ii = -1*read_k+3 438 439 call waveformat(cg,cg_disk,cg_index,phi,dk,ii,ikpt,& 440 & ikpt_,isgn,isppol,jj,jkpt,jkpt_,kg_kpt,kpt,kg_jl,maxband,mband,mcg,mcg_disk,& 441 & minband,mkmem,mpw,nkpt,nkpt_,npwarr,nsppol,nspinor,shift_g_2,tr) 442 443 ! Compute the overlap matrix <u_k|u_k+b> 444 445 do iband=minband,maxband 446 cg_index_iband=cg_index(iband,ikpt,isppol) 447 do jband=minband,maxband 448 do ipw=1,npwarr(ikpt) 449 cmatrix(1,iband,jband)=cmatrix(1,iband,jband)+& 450 & cg(1,ipw+cg_index_iband)*phi(1,ipw,jband)+& 451 & tr(ii)*cg(2,ipw+cg_index_iband)*tr(jj)*phi(2,ipw,jband) 452 453 cmatrix(2,iband,jband)=cmatrix(2,iband,jband)+& 454 & cg(1,ipw+cg_index_iband)*tr(jj)*phi(2,ipw,jband)-& 455 & tr(ii)*cg(2,ipw+cg_index_iband)*phi(1,ipw,jband) 456 end do 457 end do 458 end do 459 460 ! Compute the inverse of cmatrix(1:2,minband:maxband, minband:maxband) 461 462 band_in = maxband - minband + 1 463 ABI_ALLOCATE(ipvt,(maxband)) 464 ABI_ALLOCATE(zgwork,(2,1:maxband)) 465 466 ! Last argument of zgedi means calculate inverse only 467 call dzgefa(cmatrix(1,minband,minband),maxband, band_in,ipvt,info) 468 call dzgedi(cmatrix(1,minband,minband),maxband, band_in,ipvt,det,zgwork,01) 469 470 ABI_DEALLOCATE(zgwork) 471 ABI_DEALLOCATE(ipvt) 472 473 ! Compute the product of Inverse overlap matrix with the wavefunction 474 475 do iband=minband,maxband 476 do ipw=1,npwarr(ikpt) 477 u_tilde(1,ipw,iband,ifor)= & 478 & dot_product(cmatrix(1,minband:maxband,iband),& 479 & phi(1,ipw,minband:maxband))-& 480 & dot_product(cmatrix(2,minband:maxband,iband),& 481 & tr(jj)*phi(2,ipw,minband:maxband)) 482 u_tilde(2,ipw,iband,ifor)= & 483 & dot_product(cmatrix(1,minband:maxband,iband),& 484 & tr(jj)*phi(2,ipw,minband:maxband))+& 485 & dot_product(cmatrix(2,minband:maxband,iband),& 486 & phi(1,ipw,minband:maxband)) 487 end do 488 end do 489 ABI_DEALLOCATE(cmatrix) 490 ABI_DEALLOCATE(phi) 491 492 end do !ifor 493 494 ! Compute dudk for ikpt 495 496 npw_k=npwarr(ikpt) 497 498 do iband=minband,maxband 499 500 icg=(iband-minband)*npw_k 501 502 dudk(1,1+icg:npw_k+icg)=(u_tilde(1,1:npw_k,iband,1)-& 503 & u_tilde(1,1:npw_k,iband,2))/twodk 504 505 dudk(2,1+icg:npw_k+icg)=(u_tilde(2,1:npw_k,iband,1)-& 506 & u_tilde(2,1:npw_k,iband,2))/twodk 507 508 end do 509 510 tim_rwwf=0 511 mcg_disk=mpw*nspinor*mband 512 nband_diff=maxband-minband+1 513 call rwwf(dudk,eig_dum_2,formeig,0,0,ikpt,isppol,kg_kpt(:,:,ikpt),& 514 & mband,mcg_disk,mpi_enreg,nband_diff,nband_diff,& 515 & npw_k,nspinor,occ_dum_2,2,1,tim_rwwf,wffddk) 516 517 !call wfk_read_band_block(wfk, band_block, ikpt, isppol, sc_mode, 518 ! kg_k=kg_kpt(:,:,ikpt), cg_k=dudk, eig_k=eig_dum, occ_k=occ_dum) 519 520 ABI_DEALLOCATE(u_tilde) 521 522 end do !ikpt 523 end do !isppol 524 525 ABI_DEALLOCATE(eig_dum_2) 526 ABI_DEALLOCATE(occ_dum_2) 527 ABI_DEALLOCATE(dudk) 528 529 call WffClose(wffddk,ierr) 530 !call wfk_close(wfk) 531 532 ABI_DEALLOCATE(kg_kpt) 533 ABI_DEALLOCATE(cg_index) 534 ABI_DEALLOCATE(ikpt_dk) 535 536 end do ! iberry 537 538 ABI_DEALLOCATE(shift_g_2) 539 ABI_DEALLOCATE(kpt) 540 541 write(std_out,*) 'uderiv: exit ' 542 543 end subroutine uderiv