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