TABLE OF CONTENTS
ABINIT/berryphase [ Functions ]
NAME
berryphase
FUNCTION
This routine is called in scfcv.f to compute the electronic Berry Phase polarization and the ionic contribution to the polarization Work for nsppol=1 or 2 ,but only accept nspinor=1, and mkmem=nkpt or 0, kptopt = 2 or 3
INPUTS
atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f) bdberry(4)=band limits for Berry phase contributions, 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 istwfk(nkpt_)=input option parameter that describes the storage of wfs kberry(3,nberry)= different delta k for Berry phases, 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_ sampels half the BZ if time-reversal symetrie is used kptopt=2 when time-reversal symetrie is used =3 when time-reversal symetrie is not 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. mpw=maximum dimensioned size of npw natom=number of atoms in cell nattyp(ntypat)= # atoms of each type. nband(nkpt*nsppol)=number of bands at each k point, for each polarization nberry=number of Berry phases to be computed nkpt=number of k points npwarr(nkpt)=number of planewaves in basis at this k point nspinor=number of spinorial components (on current proc) nsppol=1 for unpolarized, 2 for spin-polarized ntypat=number of types of atoms in unit cell nkpt_=number of k points generated by ABINIT, (see kpt_) rprimd(3,3)=dimensional real space primitive translations (bohr) ucvol=unit cell volume in bohr**3. xred(3,natom)=reduced atomic coordinates zion(ntypat)=valence charge of each type of atom
OUTPUT
(the polarization is printed)
SIDE EFFECTS
TODO
Cleaning, checking for rules. Should allow for time-reversal symmetry (istwfk) Should use randac to scan rapidly the wf file
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 determinant(:,:)= determinant of cmatrix det_average(2)= averaged det_string over all strings det_string(:,:)= determinant product of cmatrices along each string dk(3)= step taken to the next k mesh point along the kberry direction dkptnext(3)= step between the next and current k point dphase= phase angle computed from rel_string(2) 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 reversal symetrie kpt_mark(nkpt)= 0, if k point is unmarked; 1, if k point has been marked maxband/minband= control the minimum and maximum band calculated in the overlap matrix nkstr= number of k points per string npw_k= npwarr(ikpt), number of planewaves in basis at this k point nstr= number of k point strings nkpt=number of k points in the whole BZ phase0= phase angle computed from det_average(2) polberry(:)= berry phase of each string (2/nsppol)*(phase0+dphase)/two_pi polb(isppol) = total berry phase polarization for each spin polbtot= total berry phase polarization polion= ionic polarization for each ion politot= total ionic polarization poltot= total polarization = polbtot + politot rel_string(2)= det_string(2)/det_average(2) shift_g(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 xcart(3,natom)= cartesian coordinates of atoms (bohr) xcart_reindex(:,:,:)= unpack xcart for each atomic species and number of atoms for each species WARNING This routine is not yet memory optimized It might be also rather time-consuming, since there is a double loop on the number of plane waves.
SOURCE
150 subroutine berryphase(atindx1,bdberry,cg,gprimd,istwfk,kberry,kg,kpt_,& 151 & kptopt,kptrlatt,mband,mcg,& 152 & mkmem,mpw,natom,nattyp,nband,nberry,npwarr,nspinor,nsppol,ntypat,& 153 & nkpt_,rprimd,ucvol,xred,zion) 154 155 !Arguments ------------------------------------ 156 !scalars 157 integer,intent(in) :: kptopt,mband,mcg,mkmem,mpw,natom,nberry,nkpt_ 158 integer,intent(in) :: nspinor,nsppol,ntypat 159 real(dp),intent(in) :: ucvol 160 !arrays 161 integer,intent(in) :: atindx1(natom),bdberry(4),istwfk(nkpt_),kberry(3,nberry) 162 integer,intent(in) :: kg(3,mpw*mkmem),kptrlatt(3,3),nattyp(ntypat) 163 integer,intent(in) :: nband(nkpt_*nsppol),npwarr(nkpt_) 164 real(dp),intent(in) :: cg(2,mcg),gprimd(1:3,1:3) 165 real(dp),intent(in) :: kpt_(3,nkpt_),rprimd(3,3),xred(3,natom),zion(ntypat) 166 167 !Local variables ------------------------- 168 !scalars 169 integer :: band_in,cg_index_iband,cg_index_jband,flag1,iatom 170 integer :: iattyp,iband,iberry,icg,ii,ikpt,ikpt2,index,index1,info 171 integer :: ipw,isppol,istr,itypat,iunmark,jband,jj,jkpt,jkstr 172 integer :: jkstr_ori,jpw,lkstr,lkstr_ori,lkstr_ori_,maxband 173 integer :: minband,nband_k,nkpt,nkstr,npw_k,nstr,read_k 174 real(dp) :: det_mod,dphase,fac,gmod,phase0,pol,polbtot,polion,politot 175 real(dp) :: poltot 176 character(len=500) :: message 177 !arrays 178 integer :: dg(3),kpt_flag(2*nkpt_),kpt_mark(2*nkpt_) 179 integer,allocatable :: cg_index(:,:,:),ikpt_dk(:),ikstr(:,:),ipvt(:) 180 integer,allocatable :: kg_dum(:,:),kg_kpt(:,:,:) 181 real(dp) :: det(2,2),det_average(2),diffk(3),dk(3),gpard(3) 182 real(dp) :: klattice(3,3),kptrlattr(3,3),polb(nsppol),rel_string(2),tr(2) 183 real(dp) :: xcart(3,natom) 184 real(dp),allocatable :: cmatrix(:,:,:),det_string(:,:) 185 real(dp),allocatable :: det_tmp(:,:),determinant(:,:),kpt(:,:) 186 real(dp),allocatable :: polberry(:),xcart_reindex(:,:,:) 187 real(dp),allocatable :: zgwork(:,:) 188 logical,allocatable :: shift_g(:) 189 190 ! *********************************************************************** 191 192 !DEBUG 193 !write(std_out,*)' berryphase : enter ' 194 !ENDDEBUG 195 196 if(nspinor==2)then 197 message = ' berryphase : does not yet work for nspinor=2' 198 ABI_ERROR(message) 199 end if 200 201 if(maxval(istwfk(:))/=1)then 202 write(message, '(a,a,a)' )& 203 & ' Sorry, this routine does not work yet with istwfk/=1.',ch10,& 204 & ' This should have been tested previously ...' 205 ABI_BUG(message) 206 end if 207 208 !change8: set up the whole k point grid in the case where kptopt = 2 209 if (kptopt==3) then 210 nkpt = nkpt_ 211 ABI_MALLOC(kpt,(3,nkpt)) 212 kpt(:,:)=kpt_(:,:) 213 else if (kptopt==2) then 214 nkpt = nkpt_*2 215 ABI_MALLOC(kpt,(3,nkpt)) 216 do ikpt = 1,nkpt/2 217 kpt_flag(ikpt) = 0 218 kpt(:,ikpt)=kpt_(:,ikpt) 219 end do 220 index = 0 221 do ikpt = (nkpt/2+1),nkpt 222 flag1 = 0 223 do jkpt = 1, nkpt/2 224 if (((abs(kpt_(1,ikpt-nkpt/2)+kpt_(1,jkpt))<1.0d-8).or.& 225 & (abs(1-abs(kpt_(1,ikpt-nkpt/2)+kpt_(1,jkpt)))<1.0d-8))& 226 & .and.((abs(kpt_(2,ikpt-nkpt/2)+kpt_(2,jkpt))<1.0d-8).or.& 227 & (abs(1-abs(kpt_(2,ikpt-nkpt/2)+kpt_(2,jkpt)))<1.0d-8))& 228 & .and.((abs(kpt_(3,ikpt-nkpt/2)+kpt_(3,jkpt))<1.0d-8).or.& 229 & (abs(1-abs(kpt_(3,ikpt-nkpt/2)+kpt_(3,jkpt)))<1.0d-8))) then 230 flag1 = 1 231 index = index + 1 232 exit 233 end if 234 end do 235 if (flag1==0) then 236 kpt_flag(ikpt-index)=ikpt-nkpt/2 237 kpt(:,ikpt-index)=-kpt_(:,ikpt-nkpt/2) 238 end if 239 end do 240 nkpt = nkpt - index 241 end if 242 243 !change8 244 245 ABI_MALLOC(shift_g,(nkpt)) 246 ABI_MALLOC(kg_dum,(3,0)) 247 248 !Compute primitive vectors of the k point lattice 249 !Copy to real(dp) 250 kptrlattr(:,:)=kptrlatt(:,:) 251 !Go to reciprocal space (in reduced coordinates) 252 call matr3inv(kptrlattr,klattice) 253 254 do iberry=1,nberry 255 256 ! Calculate dimensional recip lattice vector along which P is calculated 257 ! dk = step to the nearest k point along that direction 258 ! in reduced coordinates 259 dk(:)=kberry(1,iberry)*klattice(:,1)+& 260 & kberry(2,iberry)*klattice(:,2)+& 261 & kberry(3,iberry)*klattice(:,3) 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 kberry direction 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_MALLOC(ikpt_dk,(nkpt)) 272 shift_g(:)= .false. 273 do ikpt=1,nkpt 274 do ikpt2=1,nkpt 275 diffk(:)=abs(kpt(:,ikpt2)-kpt(:,ikpt)-dk(:)) 276 if(sum(abs(diffk(:)-nint(diffk(:))))<3*tol8)then 277 ikpt_dk(ikpt)=ikpt2 278 if(sum(diffk(:))>=3*tol8)shift_g(ikpt2) = .true. 279 exit 280 end if 281 end do 282 end do 283 284 ! DEBUG 285 ! do ikpt = 1,nkpt 286 ! write(100,*)'ikpt_dk = ',ikpt_dk(ikpt) 287 ! if (shift_g(ikpt))then 288 ! write(100,*)'true' 289 ! else 290 ! write(100,*)'false' 291 ! end if 292 ! write(100,*)'' 293 ! end do 294 ! ENDDEBUG 295 296 ! Find the string length, starting from k point 1 297 ! (all strings must have the same number of points) 298 nkstr=1 299 ikpt2=1 300 do ikpt=1,nkpt 301 ikpt2=ikpt_dk(ikpt2) 302 if(ikpt2==1)exit 303 nkstr=nkstr+1 304 end do 305 306 ! Check that the string length is a divisor of nkpt 307 if(mod(nkpt,nkstr)/=0)then 308 write(message,'(a,a,a,a,i5,a,i7)')ch10,& 309 & ' berryphase: BUG -',ch10,& 310 & ' The string length=',nkstr,', is not a divisor of nkpt=',nkpt 311 call wrtout(std_out,message,'COLL') 312 end if 313 nstr=nkpt/nkstr 314 315 write(message,'(a,a,a,3f9.5,a,a,3f9.5,a)')ch10,& 316 & ' Computing the polarization (Berry phase) for reciprocal vector:',ch10,& 317 & dk(:),' (in reduced coordinates)',ch10,& 318 & gpard(1:3),' (in cartesian coordinates - atomic units)' 319 call wrtout(ab_out,message,'COLL') 320 call wrtout(std_out,message,'COLL') 321 322 write(message,'(a,i5,a,a,i5)')& 323 & ' Number of strings: ',nstr,ch10,& 324 & ' Number of k points in string:', nkstr 325 call wrtout(std_out,message,'COLL') 326 327 if(nsppol==1)then 328 write(message, '(a,i5,a,i5)')& 329 & ' From band number',bdberry(1),' to band number',bdberry(2) 330 else 331 write(message, '(a,i5,a,i5,a,a,a,i5,a,i5,a)')& 332 & ' From band number',bdberry(1),' to band number',bdberry(2),' for spin up,',& 333 & ch10,& 334 & ' from band number',bdberry(3),' to band number',bdberry(4),' for spin down.' 335 end if 336 call wrtout(ab_out,message,'COLL') 337 call wrtout(std_out,message,'COLL') 338 339 ! DEBUG 340 ! write(std_out,*)' berryphase : find nkpt,nkstr,nstr=',nkpt,nkstr,nstr 341 ! stop 342 ! ENDDEBUG 343 344 ! Build the different strings 345 ABI_MALLOC(ikstr,(nkstr,nstr)) 346 347 iunmark=1 348 kpt_mark(:)=0 349 do istr = 1, nstr 350 do while(kpt_mark(iunmark)/=0) 351 iunmark = iunmark + 1 352 end do 353 ikstr(1, istr) = iunmark 354 kpt_mark(iunmark)=1 355 do jkstr = 2, nkstr 356 ikstr(jkstr,istr)=ikpt_dk(ikstr(jkstr-1,istr)) 357 kpt_mark(ikstr(jkstr,istr))=1 358 end do 359 end do ! istr 360 361 ! DEBUG 362 ! do istr = 1,nstr 363 ! do jkstr = 1,nkstr 364 ! if (shift_g(ikstr(jkstr,istr))) then 365 ! write(99,*) ikstr(jkstr,istr),'true' 366 ! else 367 ! write(99,*) ikstr(jkstr,istr),'false' 368 ! end if 369 ! end do 370 ! end do 371 ! ENDDEBUG 372 373 ABI_FREE(ikpt_dk) 374 ! DEBUG! 375 ! write(100,*) 'list all the k points strings:' 376 ! do istr=1,nstr 377 ! write(100,*) (ikstr(jkstr,istr),jkstr=1,nkstr) 378 ! end do 379 ! ENDDEBUG! 380 381 ! ***************************************************************************** 382 ! Find the location of each wavefunction 383 ABI_MALLOC(cg_index,(mband,nkpt,nsppol)) 384 385 icg = 0 386 do isppol=1,nsppol 387 do ikpt=1,nkpt_ 388 nband_k=nband(ikpt+(isppol-1)*nkpt_) 389 npw_k=npwarr(ikpt) 390 do iband=1,nband_k 391 cg_index(iband,ikpt,isppol)=(iband-1)*npw_k*nspinor+icg 392 end do 393 icg=icg+npw_k*nspinor*nband(ikpt) 394 end do 395 end do 396 397 ! change5 398 if (mkmem/=0) then 399 ! Find the planewave vectors and their indexes for each k point 400 ABI_MALLOC(kg_kpt,(3,mpw*nspinor,nkpt_)) 401 kg_kpt(:,:,:) = 0 402 index1 = 0 403 do ikpt=1,nkpt_ 404 npw_k=npwarr(ikpt) 405 do ipw=1,npw_k*nspinor 406 kg_kpt(1:3,ipw,ikpt)=kg(1:3,ipw+index1) 407 end do 408 index1=index1+npw_k*nspinor 409 end do 410 end if !change5 411 ! ***************************************************************************** 412 ABI_MALLOC(det_string,(2, nstr)) 413 ABI_MALLOC(det_tmp,(2, nstr)) 414 ABI_MALLOC(polberry,(nstr)) 415 416 ! Initialize berry phase polarization for each spin and the total one 417 polb(1:nsppol) = 0.0_dp 418 polbtot=0.0_dp 419 420 ! Loop over spins 421 do isppol=1,nsppol 422 423 minband=bdberry(2*isppol-1) 424 maxband=bdberry(2*isppol) 425 426 if(minband<1)then 427 write(message,'(a,i0,a)')' The band limit minband=',minband,', is lower than 0.' 428 ABI_BUG(message) 429 end if 430 431 if(maxband<1)then 432 write(message,'(a,i0,a)')' The band limit maxband=',maxband,', is lower than 0.' 433 ABI_BUG(message) 434 end if 435 436 if(maxband<minband)then 437 write(message,'(a,i0,a,i0)')' maxband=',maxband,', is lower than minband=',minband 438 ABI_BUG(message) 439 end if 440 441 ! Initialize det_string and det_average 442 det_string(1, 1:nstr) = 1.0_dp; det_string(2, 1:nstr) = 0.0_dp 443 det_average(1:2)=0.0_dp; det_average(2)=0.0_dp 444 445 ! Loop over strings 446 do istr = 1, nstr 447 448 ! change7 449 read_k = 0 450 451 ! DEBUG! 452 ! write(100,'(a,i4)') 'This is in string', istr 453 ! ENDDEBUG! 454 455 ! Loop over k points per string 456 ABI_MALLOC(determinant,(2, nkstr)) 457 458 do jkstr=1,nkstr 459 460 ABI_MALLOC(cmatrix,(2,maxband,maxband)) 461 if(jkstr < nkstr) then 462 lkstr=jkstr+1 463 else 464 lkstr= jkstr+1-nkstr 465 end if 466 jkstr_ori=ikstr(jkstr,istr) 467 lkstr_ori=ikstr(lkstr,istr) 468 469 ! change9 470 lkstr_ori_=lkstr_ori 471 tr(1) = 1.0_dp 472 tr(2) = 1.0_dp 473 if (kptopt==2) then 474 if (read_k == 0) then 475 if (kpt_flag(jkstr_ori)/=0) then 476 tr(1) = -1.0_dp 477 jkstr_ori = kpt_flag(jkstr_ori) 478 end if 479 if (kpt_flag(lkstr_ori)/=0) then 480 tr(2) = -1.0_dp 481 lkstr_ori = kpt_flag(lkstr_ori) 482 end if 483 else !read_k 484 if (kpt_flag(jkstr_ori)/=0) then 485 tr(-1*read_k+3) = -1.0_dp 486 jkstr_ori = kpt_flag(jkstr_ori) 487 end if 488 if (kpt_flag(lkstr_ori)/=0) then 489 tr(read_k) = -1.0_dp 490 lkstr_ori = kpt_flag(lkstr_ori) 491 end if 492 end if !read_k 493 end if !kptopt 494 ! change9 495 496 nband_k=nband(jkstr_ori+(isppol-1)*nkpt_) 497 if(nband_k<maxband)then 498 write(message,'(a,i0,a,i0)')' maxband=',maxband,', is larger than nband(j,isppol)=',nband_k 499 ABI_BUG(message) 500 end if 501 502 nband_k=nband(lkstr_ori+(isppol-1)*nkpt_) 503 if(nband_k<maxband)then 504 write(message,'(a,i0,a,i0)')& 505 & ' maxband=',maxband,', is larger than nband(l,isppol)=',nband_k 506 ABI_BUG(message) 507 end if 508 509 if (jkstr==1) read_k = 2 510 ! Compute the overlap matrix <u_k|u_k+b> 511 cmatrix(1:2,1:maxband,1:maxband)=zero 512 jj = read_k 513 ii = -1*read_k+3 514 if(.not. shift_g(lkstr_ori_) ) then 515 ! Change3 516 do ipw=1,npwarr(jkstr_ori) 517 do jpw=1,npwarr(lkstr_ori) 518 519 ! Check if Fourier components of jkstr and jkstr+1 matches 520 521 if((tr(ii)*kg_kpt(1,ipw,jkstr_ori)==tr(jj)*kg_kpt(1,jpw,lkstr_ori))& 522 & .and.(tr(ii)*kg_kpt(2,ipw,jkstr_ori) == tr(jj)*kg_kpt(2,jpw,lkstr_ori))& 523 & .and.(tr(ii)*kg_kpt(3,ipw,jkstr_ori) == tr(jj)*kg_kpt(3,jpw,lkstr_ori)))& 524 & then 525 526 do iband=minband,maxband 527 cg_index_iband=cg_index(iband,jkstr_ori,isppol) 528 do jband=minband,maxband 529 cg_index_jband=cg_index(jband,lkstr_ori,isppol) 530 531 cmatrix(1,iband,jband)=cmatrix(1,iband,jband)+& 532 & cg(1,ipw+cg_index_iband)*cg(1,jpw+cg_index_jband)+& 533 & tr(ii)*cg(2,ipw+cg_index_iband)*tr(jj)*cg(2,jpw+cg_index_jband) 534 cmatrix(2,iband,jband)=cmatrix(2,iband,jband)+& 535 & cg(1,ipw+cg_index_iband)*tr(jj)*cg(2,jpw+cg_index_jband)-& 536 & tr(ii)*cg(2,ipw+cg_index_iband)*cg(1,jpw+cg_index_jband) 537 538 end do !jband 539 end do !iband 540 exit !stop loop over jpw if Fourier components of jkstr and jkstr + 1 matches 541 end if 542 543 end do ! jpw 544 end do ! ipw 545 546 ! But there is a special pair of k points which involves the shift of a 547 ! G vector 548 549 else 550 551 dg(:) = -1*nint(tr(jj)*kpt(:,lkstr_ori)-tr(ii)*kpt(:,jkstr_ori)-dk(:)) 552 553 ! DEBUG 554 ! write(100,*)dg 555 ! write(100,*)kberry(:,iberry) 556 ! write(100,*)'' 557 ! ENDDEBUG 558 559 ! change4 560 do ipw=1,npwarr(jkstr_ori) 561 do jpw=1,npwarr(lkstr_ori) 562 563 ! Check if Fourier components of jkstr and jkstr+1 564 ! matches by comparing the G vectors 565 566 if((tr(ii)*kg_kpt(1,ipw,jkstr_ori)==tr(jj)*kg_kpt(1,jpw,lkstr_ori)-dg(1))& 567 & .and.(tr(ii)*kg_kpt(2,ipw,jkstr_ori) == tr(jj)*kg_kpt(2,jpw,lkstr_ori)-dg(2))& 568 & .and.(tr(ii)*kg_kpt(3,ipw,jkstr_ori) == tr(jj)*kg_kpt(3,jpw,lkstr_ori)-dg(3)))& 569 & then 570 571 do iband=minband,maxband 572 cg_index_iband=cg_index(iband,jkstr_ori,isppol) 573 574 do jband=minband,maxband 575 cg_index_jband=cg_index(jband,lkstr_ori,isppol) 576 577 cmatrix(1,iband,jband)=cmatrix(1,iband,jband)+& 578 & cg(1,ipw+cg_index_iband)*cg(1,jpw+cg_index_jband)+& 579 & tr(ii)*cg(2,ipw+cg_index_iband)*tr(jj)*cg(2,jpw+cg_index_jband) 580 cmatrix(2,iband,jband)=cmatrix(2,iband,jband)+& 581 & cg(1,ipw+cg_index_iband)*tr(jj)*cg(2,jpw+cg_index_jband)-& 582 & tr(ii)*cg(2,ipw+cg_index_iband)*cg(1,jpw+cg_index_jband) 583 584 end do ! jband 585 end do ! iband 586 exit !stop loop over jpw if Fourier components of jkstr and jkstr + 1 matches 587 end if 588 end do ! jpw 589 end do ! ipw 590 end if 591 592 ! Compute the determinant of cmatrix(1:2,minband:maxband, minband:maxband) 593 594 band_in = maxband - minband + 1 595 596 ABI_MALLOC(ipvt,(maxband)) 597 ABI_MALLOC(zgwork,(2,1:maxband)) 598 599 ! Last argument of zgedi means calculate determinant only. 600 call dzgefa(cmatrix(1,minband,minband),maxband, band_in,ipvt,info) 601 call dzgedi(cmatrix(1,minband,minband),maxband, band_in,ipvt,det,zgwork,10) 602 603 ABI_FREE(zgwork) 604 ABI_FREE(ipvt) 605 606 fac=exp(log(10._dp)*det(1,2)) 607 determinant(1, jkstr) = fac*(det(1,1)*cos(log(10._dp)*det(2,2)) - & 608 & det(2,1)*sin(log(10._dp)*det(2,2))) 609 determinant(2, jkstr) = fac*(det(1,1)*sin(log(10._dp)*det(2,2)) + & 610 & det(2,1)*cos(log(10._dp)*det(2,2))) 611 ! DEBUG! 612 ! write(100,*) 'det',jkstr,lkstr,'=', determinant(1:2,jkstr) 613 ! ENDDEBUG! 614 615 det_tmp(1,istr) = det_string(1,istr)*determinant(1,jkstr) - & 616 & det_string(2,istr)*determinant(2,jkstr) 617 det_tmp(2,istr) = det_string(1,istr)*determinant(2,jkstr) + & 618 & det_string(2,istr)*determinant(1,jkstr) 619 det_string(1:2,istr) = det_tmp(1:2,istr) 620 621 ABI_FREE(cmatrix) 622 623 ! Close loop over k points along string 624 read_k = -1*read_k + 3 ! read_k=2 <-> read_k=1 625 end do 626 627 ! DEBUG! 628 ! write(100,*) 'det_string =', det_string(1:2,istr) 629 ! write(100,*) 630 ! ENDDEBUG! 631 632 det_average(1) = det_average(1) + det_string(1,istr)/nstr 633 det_average(2) = det_average(2) + det_string(2,istr)/nstr 634 635 ABI_FREE(determinant) 636 637 ! Close loop over strings 638 end do 639 640 641 ! ***************************************************************************** 642 ! Calculate the electronic contribution to the polarization 643 644 write(message,'(a,a)')ch10,& 645 & ' Compute the electronic contribution to polarization' 646 call wrtout(std_out,message,'COLL') 647 648 ! First berry phase that corresponds to det_average 649 phase0 = atan2(det_average(2),det_average(1)) 650 det_mod = det_average(1)**2+det_average(2)**2 651 652 ! Then berry phase that corresponds to each string relative to the average 653 do istr = 1, nstr 654 rel_string(1) = (det_string(1,istr)*det_average(1) + & 655 det_string(2,istr)*det_average(2))/det_mod 656 rel_string(2) = (det_string(2,istr)*det_average(1) - & 657 det_string(1,istr)*det_average(2))/det_mod 658 dphase = atan2(rel_string(2),rel_string(1)) 659 polberry(istr) = (2.0_dp/nsppol)*(phase0+dphase)/two_pi 660 polb(isppol) = polb(isppol) + polberry(istr)/nstr 661 end do 662 663 ! Output berry phase polarization 664 write(message,'(a,10x,a,10x,a)')ch10,& 665 & 'istr','polberry(istr)' 666 call wrtout(std_out,message,'COLL') 667 do istr=1,nstr 668 write(message,'(10x,i4,7x,e16.9)')istr,polberry(istr) 669 call wrtout(std_out,message,'COLL') 670 end do 671 672 write(message,'(9x,a,7x,e16.9,1x,a,i4,a,a)')& 673 & 'total',polb(isppol),'(isppol=',isppol,')',ch10 674 call wrtout(std_out,message,'COLL') 675 676 polbtot=polbtot+polb(isppol) 677 678 end do ! isppol 679 680 ABI_FREE(polberry) 681 ABI_FREE(det_tmp) 682 ABI_FREE(det_string) 683 ABI_FREE(ikstr) 684 ABI_FREE(cg_index) 685 ! change6 686 if (mkmem /=0) then 687 ABI_FREE(kg_kpt) 688 end if 689 ! ***************************************************************************** 690 ! Reindex xcart according to atom and type 691 call xred2xcart(natom,rprimd,xcart,xred) 692 ABI_MALLOC(xcart_reindex,(3,natom,ntypat)) 693 index=1 694 do itypat=1,ntypat 695 do iattyp=1,nattyp(itypat) 696 iatom=atindx1(index) 697 xcart_reindex(1:3,iattyp,itypat) = xcart(1:3,iatom) 698 index = index+1 699 end do 700 end do 701 702 ! Compute the ionic contribution to the polarization 703 politot = 0.0_dp 704 write(message,'(a)')' Compute the ionic contributions' 705 call wrtout(std_out,message,'COLL') 706 707 write(message,'(a,2x,a,2x,a,15x,a)')ch10,& 708 & 'itypat', 'iattyp', 'polion' 709 call wrtout(std_out,message,'COLL') 710 711 do itypat=1,ntypat 712 do iattyp=1,nattyp(itypat) 713 polion=zion(itypat)*nkstr*& 714 & dot_product(xcart_reindex(1:3,iattyp,itypat),gpard(1:3)) 715 ! Fold into interval (-1,1) 716 polion=polion-2._dp*nint(polion/2.0_dp) 717 politot=politot+polion 718 write(message,'(2x,i2,5x,i2,10x,e16.9)') itypat,iattyp,polion 719 call wrtout(std_out,message,'COLL') 720 end do 721 end do 722 723 ! Fold into interval (-1,1) again 724 politot=politot-2.0_dp*nint(politot/2.0_dp) 725 726 write(message,'(9x,a,7x,es19.9)') 'total',politot 727 call wrtout(std_out,message,'COLL') 728 729 ABI_FREE(xcart_reindex) 730 731 ! Compute the total polarizations 732 733 poltot=politot+polbtot 734 735 write(message,'(a,a)')ch10,& 736 & ' Summary of the results' 737 call wrtout(std_out,message,'COLL') 738 call wrtout(ab_out,message,'COLL') 739 740 write(message,'(a,es19.9)')& 741 & ' Electronic Berry phase ' ,polbtot 742 call wrtout(std_out,message,'COLL') 743 call wrtout(ab_out,message,'COLL') 744 745 write(message,'(a,es19.9)') & 746 & ' Ionic phase ', politot 747 call wrtout(std_out,message,'COLL') 748 call wrtout(ab_out,message,'COLL') 749 750 write(message,'(a,es19.9)') & 751 & ' Total phase ', poltot 752 call wrtout(std_out,message,'COLL') 753 call wrtout(ab_out,message,'COLL') 754 755 poltot=poltot-2.0_dp*nint(poltot/2._dp) 756 write(message,'(a,es19.9)') & 757 & ' Remapping in [-1,1] ', poltot 758 call wrtout(std_out,message,'COLL') 759 call wrtout(ab_out,message,'COLL') 760 761 ! Transform the phase into a polarization 762 fac = 1._dp/(gmod*nkstr) 763 fac = fac/ucvol 764 pol = fac*poltot 765 766 write(message,'(a,a,es19.9,a,a,a,es19.9,a,a)')ch10,& 767 & ' Polarization ', pol,' (a.u. of charge)/bohr^2',ch10,& 768 & ' Polarization ', pol*(e_Cb)/(Bohr_Ang*1d-10)**2,& 769 & ' C/m^2',ch10 770 call wrtout(std_out,message,'COLL') 771 call wrtout(ab_out,message,'COLL') 772 773 end do ! iberry 774 775 ABI_FREE(shift_g) 776 ABI_FREE(kpt) 777 ABI_FREE(kg_dum) 778 779 end subroutine berryphase
ABINIT/m_berryphase [ Modules ]
NAME
m_berryphase
FUNCTION
COPYRIGHT
Copyright (C) 2000-2022 ABINIT group (NSAI,XG,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
15 #if defined HAVE_CONFIG_H 16 #include "config.h" 17 #endif 18 19 #include "abi_common.h" 20 21 module m_berryphase 22 23 use defs_basis 24 use m_errors 25 use m_abicore 26 use m_hdr 27 use m_dtset 28 29 use m_geometry, only : xred2xcart 30 use m_hide_lapack, only : dzgedi, dzgefa 31 use m_symtk, only : matr3inv 32 33 implicit none 34 35 private