TABLE OF CONTENTS


ABINIT/berryphase [ Functions ]

[ Top ] [ 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.

PARENTS

      elpolariz

CHILDREN

      dzgedi,dzgefa,matr3inv,wrtout,xred2xcart

SOURCE

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

ABINIT/m_berryphase [ Modules ]

[ Top ] [ Modules ]

NAME

  m_berryphase

FUNCTION

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 .

PARENTS

CHILDREN

SOURCE

19 #if defined HAVE_CONFIG_H
20 #include "config.h"
21 #endif
22 
23 #include "abi_common.h"
24 
25 module m_berryphase
26 
27  use defs_basis
28  use defs_abitypes
29  use m_errors
30  use m_abicore
31  use m_hdr
32 
33  use m_geometry,     only : xred2xcart
34  use m_hide_lapack,  only : dzgedi, dzgefa
35  use m_symtk,        only : matr3inv
36 
37  implicit none
38 
39  private