[Wannier] guiding centers

[Nicola Marzari]

On 15/11/2013 12:07, Elias Assmann wrote:
> Hi List,
>
> I would like to understand how the “guiding centers“ that can be
> specified in the ‘win’ file work.  I have checked the user guide,
> tutorial, and review article, but cannot find anything definite.
>
> 	Elias
>


Dear Elias,

the idea is that if M_nn_kb was dispersionless, it would be identical
to b scalar rcenter, and (see eq 28 of the 2012 review) would give you 
rcenter. Since once the functions are maximally localized,
M_nn is optimally smooth as you surf the Brillouin zone, you want
to make sure that as you move around, the branch cuts are such that
M_nn gives you a Im ln M_nn that is roughly always b scalar rguide,
where rguide is an educated guess to rcenter.

That info is still buried in the wann_phases subroutine inside 
wannierise.f90 - a bit more info in this old f77 version.

			nicola



-- 

----------------------------------------------------------------------
Prof Nicola Marzari, Chair of Theory and Simulation of Materials, EPFL
-------------- next part --------------
c-----------------------------------------------------------------------

      subroutine phases(nbans,nkpts2,nntot,nnh,neigh,
     &           bk,bka,cm,csheet,sheet,rguide,irguide)
      
      INCLUDE 'wannier.include'
      IMPLICIT COMPLEX (C)
      dimension cm(nbands,nbands,nkpts,nnmx), csum(nnmxh)
      dimension bk(3,nkpts,nnmx), bka(3,nnmxh), nntot(nkpts)
      dimension neigh(nkpts,nnmxh), rguide(3,nbands)
      dimension xx(nnmx), smat(3,3),svec(3),sinv(3,3)
      dimension csheet(nbands,nkpts,nnmx),sheet(nbands,nkpts,nnmx)

      dimension rnkb(nbands,nkpts,nnmx)

      ci=(0.,1.)

c report problem to solve
c for each band, csum is determined and then its appropriate
c guiding center rguide(3,nwann)

      do nwann=1,nbans

c get average phase for each unique bk direction

       do na=1,nnh
        csum(na)=(0.0,0.0)
        do nkp=1,nkpts2
         nn=neigh(nkp,na)
         csum(na)=csum(na)+cm(nwann,nwann,nkp,nn)
        end do
       end do

c now analyze that information to get good guess at
c wannier center

c      write(*,*)
c      do na=1,nnh
c       write(*,'a,3f10.5,a,2f10.5)')
c    &    ' bka=',(bka(j,na),j=1,3),'  csum=',csum(na)
c      end do

c problem is to find a real-space 3-vector rguide such that
c   phase of csum(nn) ~= phase of exp[ -i bka(nn) dot rguide ]
c or, letting
c   xx(nn) = - Im ln csum(nn)  (modulo 2*pi)
c then
c   bka(nn) dot rguide ~= xx(nn)
c
c we take an arbitrary branch cut for first three xx(nn)
c and determine rguide from these; then for each additional bka
c vector, we first determine the most consistent branch cut,
c and then update rguide
c
c rguide is obtained by minimizing
c   sum_nn [ bka(nn) dot rguide - xx(nn) ] ^2
c or, setting the derivative with respect to rcenter to zero,
c   sum_i smat(j,i) * rguide(i,nwann) = svec(j)
c where
c   smat(j,i) = sum_nn bka(j,nn) * bka(i,nn)
c   svec(j)   = sum_nn bka(j,nn) * xx(nn)

c initialize smat and svec
       do j=1,3
        do i=1,3
          smat(j,i)=0.
        end do
        svec(j)=0.
       end do

       write(*,*) ' '

       do nn=1,nnh

        if (nn.le.3) then
c         obtain xx with arbitrary branch cut choice
          xx(nn)=-dimag(log(csum(nn)))
        else
c         obtain xx with branch cut choice guided by rguide
          xx0=0.
          do j=1,3
            xx0=xx0+bka(j,nn)*rguide(j,nwann)
          end do
c         xx0 is expected value for xx
          csumt=cexp(ci*xx0)
c         csumt has opposite of expected phase of csum(nn)
          xx(nn)=xx0-dimag(log(csum(nn)*csumt))
        endif

        write(*,'(a,i5,3f7.3,2f10.5)') 'nn, bka, xx, mag =',
     1    nn,(bka(j,nn),j=1,3),xx(nn),abs(csum(nn))/float(nkpts2)

c       update smat and svec
        do j=1,3
          do i=1,3
            smat(j,i)=smat(j,i)+bka(j,nn)*bka(i,nn)
          end do
          svec(j)=svec(j)+bka(j,nn)*xx(nn)
        end do

        if (nn.ge.3) then
c         determine rguide
          call inv3(smat,sinv,det)
c         the inverse of smat is sinv/det
          if (abs(det).lt.1.e-06) then
c           evidently first three bka vectors are linearly dependent
c           this is not allowed
            write(*,*) ' error in findr: dependency'
            stop
          endif
          if (irguide.ne.0) then
          do j=1,3
            rguide(j,nwann)=0.
            do i=1,3
              rguide(j,nwann)=rguide(j,nwann)+sinv(j,i)*svec(i)/det
            end do
          end do
          end if
          write(*,'(a,3f10.5)') 'rguide =',(rguide(i,nwann),i=1,3)
        endif

       end do

      end do

c     obtain branch cut choice guided by rguid
  
      do nkp=1,nkpts2
       do nwann=1,nbans
        do nn=1,nntot(nkp)
         sheet(nwann,nkp,nn)=0.0
         do j=1,3
          sheet(nwann,nkp,nn)=sheet(nwann,nkp,nn)+
     &                        bk(j,nkp,nn)*rguide(j,nwann)
         end do
         csheet(nwann,nkp,nn)=cexp(ci*sheet(nwann,nkp,nn))
        end do
       end do
      end do

c now check that we picked the proper sheet for the log
c of cm. criterion: q_n^{k,b}=Im(ln(M_nn^{k,b})) + b \cdot r_n are
c circa 0 for a good solution, circa multiples of 2 pi  for a bad one.
c I use the guiding center, instead of r_n, to understand which could be
c right sheet

c     write(*,*) ' '
c     write(*,*) 'Approximate q_n^{k,b}, using the guiding centers'
c     write(*,*) ' '

      do nkp=1,nkpts2
       do m=1,nbans
        do nn=1,nntot(nkp)
         rnkb(m,nkp,nn)=0.0
         brn=0.0
         do ind=1,3
          brn=brn+bk(ind,nkp,nn)*rguide(ind,m)
         end do
         rnkb(m,nkp,nn)=rnkb(m,nkp,nn)+brn
        end do
       end do
      end do
      do nkp=1,nkpts2
       do n=1,nbans
        do nn=1,nntot(nkp)
         pherr=dimag(log(csheet(n,nkp,nn)*cm(n,n,nkp,nn)))
     &   -sheet(n,nkp,nn)+rnkb(n,nkp,nn)-dimag(log(cm(n,n,nkp,nn)))
         if (abs(pherr).gt.3.1416) then
      write(*,'(3i4,f18.9,3f10.5)') nkp,n,nn,pherr,(bk(i,nkp,nn),i=1,3)
         end if
        end do
100     continue
       end do
      end do

      return
      end