[Wannier] position matrix (dipole transition matrix) printed out in wannier

[Ivo Souza]

Dear Liping,

  There is one other place in the wannier90 package where the
Wannier-basis position matrix elements <0m|r|Rn> are calculated: the
subroutine get_AA_R in src/postw90/get_oper.F90

  (Those matrix elements are extensively used in the "berry" module of
postw90, because their Fourier transform gives the so-called "Berry
connection matrix" A_{mn}(k).)

  You might want to take a look into that subroutine, and see if you can
re-use it for your purposes.

  Best,

  Ivo

On Wed, 15 Jan 2014, Liping Chen wrote:

> Dear Wannier users,
>
> I need to get the info of dipole transition matrix <i|er|j>=e<i|r|j>. So I just need to get the position matrix <i|r|j>. But it seems that there is no way to ask wannier to print out the position matrix in seedname.win. I find that the matrix elements of r2 between WF can be printed out by setting WRITE_R2MN = T. And the codes for writing the matrix elements of r2 are like this:
> ###########################################################################
>     open(r2mnunit,file=trim(seedname)//'.r2mn',form='formatted',err=158)
>     do nw1 = 1, num_wann
>        do nw2 = 1, num_wann
>           r2ave_mn = 0.0_dp
>           delta = 0.0_dp
>           if (nw1.eq.nw2) delta = 1.0_dp
>           do nkp = 1, num_kpts
>              do nn = 1, nntot
>                 r2ave_mn = r2ave_mn + wb(nn) * &
>                      ! [GP-begin, Apr13, 2012: corrected sign inside "real"]
>                      ( 2.0_dp * delta - real(m_matrix(nw1,nw2,nn,nkp) + &
>                      conjg(m_matrix(nw2,nw1,nn,nkp)),kind=dp) )
>                      ! [GP-end]
>              enddo
>           enddo
>           r2ave_mn = r2ave_mn / real(num_kpts,dp)
>           write (r2mnunit, '(2i6,f20.12)') nw1, nw2, r2ave_mn
>        enddo
>     enddo
>     close(r2mnunit)
> ###########################################################################
>
> There should be some similarities for calculating r2ave and the above r2. I've found the codes to calculate rave, rave2 and r2ave in subroutine wann_omega.
> ##########################################################################################
>     do nkp = 1, num_kpts
>        do nn = 1, nntot
>           do n = 1, num_wann
>              ! Note that this ln_tmp is defined differently wrt the one in wann_domega
>              ln_tmp(n,nn,nkp)=( aimag(log(csheet(n,nn,nkp) &
>                      * m_matrix(n,n,nn,nkp))) - sheet(n,nn,nkp) )
>           end do
>       end do
>     end do
>
>
>     rave  = 0.0_dp
>     do iw = 1, num_wann
>        do ind = 1, 3
>           do nkp = 1, num_kpts
>              do nn = 1, nntot
>                 rave(ind,iw) = rave(ind,iw) + wb(nn) * bk(ind,nn,nkp) &
>                       *ln_tmp(iw,nn,nkp)
>              enddo
>           enddo
>        enddo
>     enddo
>     rave = -rave/real(num_kpts,dp)
>
>     rave2 = 0.0_dp
>     do iw = 1, num_wann
>        rave2(iw) = sum(rave(:,iw)*rave(:,iw))
>     enddo
>
>     r2ave = 0.0_dp
>     do iw = 1, num_wann
>        do nkp = 1, num_kpts
>           do nn = 1, nntot
>              mnn2 = real(m_matrix(iw,iw,nn,nkp)*conjg(m_matrix(iw,iw,nn,nkp)),kind=dp)
>              r2ave(iw) = r2ave(iw) + wb(nn) * ( 1.0_dp - mnn2 + ln_tmp(iw,nn,nkp)**2 )
>           enddo
>        enddo
>     enddo
>     r2ave = r2ave/real(num_kpts,dp)
> ##########################################################################################
> It's true the calculation of r2ave and r2 are similar. But why there is an extra term of "ln_tmp(iw,nn,nkp)**2" for r2ave but not for r2?
>
> What should I do to calculate r?
> ###########################################################################
>     open(rmnunit,file=trim(seedname)//'.rmn',form='formatted',err=158)
>     do nw1 = 1, num_wann
>        do nw2 = 1, num_wann
>           rave_mn = 0.0_dp
>           delta = 0.0_dp
>           if (nw1.eq.nw2) delta = 1.0_dp
>           do nkp = 1, num_kpts
>              do nn = 1, nntot
>                 rave_mn = rave_mn + wb(nn) * ?????? (What should I put here?)
>              enddo
>           enddo
>           rave_mn = rave_mn / real(num_kpts,dp)
>           write (rmnunit, '(2i6,f20.12)') nw1, nw2, rave_mn
>        enddo
>     enddo
>     close(rmnunit)
> ###########################################################################
>
> Thanks,
>
> Liping Chen
> Department of Chemistry
> University of Rochester
> _______________________________________________
> Wannier mailing list
> Wannier at quantum-espresso.org
> http://www.democritos.it/mailman/listinfo/wannier
>