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

[Liping Chen]

Dear Ivo,

Thank you very much for your reply. I've had a look at the subroutine get_AA_R. I want to make sure with you that near the end of the subroutine there are
############################################################
       call fourier_q_to_R(AA_q(:,:,:,1),AA_R(:,:,:,1))
       call fourier_q_to_R(AA_q(:,:,:,2),AA_R(:,:,:,2))
       call fourier_q_to_R(AA_q(:,:,:,3),AA_R(:,:,:,3))
############################################################
AA_R(:,:,:,1), AA_R(:,:,:,2) and AA_R(:,:,:,3) are the position matrix elements <0m|r|Rn> I want, right?

Here, I also have a question about the idir = 1, 2, 3. Are these three directions along the three lattice vectors? If so, what about when the lattice vectors are not orthogonal？

Thanks,

Liping 

----- Original Message -----
From: "Ivo Souza" <ivo_souza at ehu.es>
To: "Liping Chen" <liping.chen at rochester.edu>
Cc: wannier at quantum-espresso.org
Sent: Thursday, January 16, 2014 3:43:10 AM
Subject: Re: [Wannier] position matrix (dipole transition matrix) printed out in wannier

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
> https://urldefense.proofpoint.com/v1/url?u=http://www.democritos.it/mailman/listinfo/wannier&k=p4Ly7qpEBiYPBVenR9G2iQ%3D%3D%0A&r=MhxMznMChIKt4IuapbjnvwFAOKg44L9%2FDsTFLScyNxw%3D%0A&m=LrMVtMdPoy5Xhv0RjvETLxTHjF8P6qJE8zOW9sqIC2w%3D%0A&s=d94ca4d5db473173648561905709bdd4b2f1c6d9fdd616d6d14989568372a2ec
>