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.

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 ]

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