
! Copyright (C) 2015 Manh Duc Le, 2017-18 Arsenii Gerasimov, Yaroslav Kvashnin
! and Lars Nordstrom. This file is distributed under the terms of the GNU
! General Public License. See the file COPYING for license details.

subroutine writew90mmn
use modmain
use modw90
implicit none
! local variables
integer ik,jk,ist,jst,i
integer ngp(nspnfv),ngpq(nspnfv)
real(8) vl(3),vc(3),q
real(8) vkql(3),t1
character(256) fname
! automatic arrays
complex(8) ylmq(lmmaxo),sfacq(natmtot)
! allocatable arrays
integer, allocatable :: igpig(:,:),igpqig(:,:)
real(8), allocatable :: jlqr(:,:)
complex(8), allocatable :: wfmt(:,:,:,:),wfir(:,:,:)
complex(8), allocatable :: wfmtq(:,:,:,:),wfgpq(:,:,:)
complex(8), allocatable :: expqmt(:,:),oq(:,:)
! find the maximum size of the spherical Bessel function array over all species
call findnjcmax
! allocate local arrays
allocate(igpig(ngkmax,nspnfv),igpqig(ngkmax,nspnfv))
allocate(jlqr(njcmax,nspecies))
allocate(wfmt(npcmtmax,natmtot,nspinor,num_bands))
allocate(wfir(ngtot,nspinor,num_bands))
allocate(wfmtq(npcmtmax,natmtot,nspinor,num_bands))
allocate(wfgpq(ngkmax,nspinor,num_bands))
allocate(expqmt(npcmtmax,natmtot))
allocate(oq(num_bands,num_bands))
fname=trim(seedname)//'.mmn'
open(50,file=trim(fname),action='WRITE',form='FORMATTED')
write(50,'("Generated by Elk version ",I1.1,".",I1.1,".",I2.2)') version
write(50,'(3I8)') num_bands,nkptnr,nntot
t1=sqrt(omega)
do ik=1,nkptnr
  call genwfsvp(.false.,.false.,num_bands,idxw90,ngridg,igfft,vkl(:,ik),ngp, &
   igpig,wfmt,ngtot,wfir)
! multiply the interstitial wavefunction by sqrt(omega)
  wfir(:,:,:)=t1*wfir(:,:,:)
  do i=1,nntot
    jk=nnlist(ik,i)
! q-vector in lattice coordinates
    vl(:)=dble(nncell(:,ik,i))+vkl(:,jk)-vkl(:,ik)
! q-vector in Cartesian coordinates
    call r3mv(bvec,vl,vc)
! q-vector length
    q=sqrt(vc(1)**2+vc(2)**2+vc(3)**2)
! generate phase factor function exp(iq.r) in the muffin-tins
    call genjlgpr(1,q,jlqr)
    call genylmv(lmaxo,vc,ylmq)
    call gensfacgp(1,vc,1,sfacq)
    call genexpmt(1,jlqr,ylmq,1,sfacq,expqmt)
! k+q-vector in lattice coordinates
    vkql(:)=vkl(:,ik)+vl(:)
! generate the wavefunctions at k+q
    call genwfsvp(.false.,.true.,num_bands,idxw90,ngridg,igfft,vkql,ngpq, &
     igpqig,wfmtq,ngkmax,wfgpq)
! determine the overlap matrix
    call genolpq(num_bands,expqmt,ngpq,igpqig,wfmt,wfir,wfmtq,wfgpq,oq)
! write overlap matrix to file
    write(50,'(5I8)') ik,jk,nncell(:,ik,i)
    do jst=1,num_bands
      do ist=1,num_bands
        write(50,'(2G18.10)') dble(oq(jst,ist)),-aimag(oq(jst,ist))
      end do
    end do
  end do
end do
close(50)
write(*,*)
write(*,'("Info(writew90mmn): created file ",A)') trim(fname)
deallocate(igpig,igpqig,jlqr)
deallocate(wfmt,wfir,wfmtq,wfgpq)
deallocate(expqmt,oq)
return
end subroutine

