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

[Liping Chen]

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