LCOV - code coverage report
Current view: top level - wannier - wann_gwf_commat.f (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 264 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 4 0.0 %

          Line data    Source code
       1             : c*************************c
       2             : c  routine to set up the  c
       3             : c  composite matrices as  c
       4             : c  input for wannier90    c
       5             : c*************************c
       6             :       module m_wann_gwf_commat
       7             :       USE m_juDFT
       8             :       implicit none
       9             :       contains
      10             : 
      11           0 :       subroutine wann_gwf_commat(nkpts,nntot,bpt,nqpts,nntot_q,
      12           0 :      >                           bpt_q,gb,gb_q,latt_const_q,
      13             :      >                           l_unformatted,l_amn,l_mmn,l_dim,
      14           0 :      >                           nparampts,param_vec)
      15             :       use m_wann_gwf_tools
      16             :       use m_wann_gwf_auxovlp
      17             :       use m_constants, only : pimach
      18             :       use m_intgr, only : intgz0
      19             : 
      20             :       implicit none
      21             :       
      22             :       ! input parameters
      23             :       logical,intent(in) :: l_unformatted,l_amn,l_mmn,l_dim(3)
      24             :       integer,intent(in) :: nkpts,nntot,nqpts,nntot_q,nparampts
      25             :       integer,intent(in) :: bpt(nntot,nkpts),bpt_q(nntot_q,nqpts)
      26             :       integer,intent(in) :: gb(3,nntot,nkpts),gb_q(3,nntot_q,nqpts)
      27             :       real,   intent(in) :: latt_const_q,param_vec(3,nparampts)
      28             : 
      29             :       ! allocatable arrays
      30           0 :       integer,allocatable :: gb_kq(:,:,:)
      31           0 :       integer,allocatable :: bpt_kq(:,:)
      32           0 :       integer,allocatable :: g1(:,:,:),g2(:,:)
      33           0 :       integer,allocatable :: b1(:,:),b2(:)
      34           0 :       real, allocatable :: xxr(:),xxi(:)
      35           0 :       complex,allocatable :: tmp_a(:,:,:),tmp_m(:,:,:,:)
      36           0 :       complex,allocatable :: mmnk(:,:,:,:),mmnq(:,:,:,:)
      37           0 :       complex,allocatable :: mk(:,:,:),mq(:,:,:)
      38           0 :       complex,allocatable :: mmn(:,:,:,:)
      39           0 :       complex,allocatable :: amn_arti2(:)
      40           0 :       character(len=20),allocatable :: feig(:),famn(:)
      41           0 :       character(len=20),allocatable :: fmmn(:),fmmn2(:)
      42             : 
      43             :       ! further variables
      44             :       integer :: nbnd,nwfs,kb,qb,nn_kq,kqb,arr_len,shift(3)
      45             :       integer :: i,j,ikqpt,ikqpt_b,nk,nq,g(3)
      46             :       integer :: nbnd_arti,nqpts_arti,nntot_arti
      47             :       integer :: nkqpts,nntot_kq,nn,ikqpt_help
      48             :       integer :: dump,nn_arti,nkp,k,b,q,nwf,n,is,ie,nx
      49             :       real :: mmn_r,mmn_i,tau,tpi,a_arti,b_arti,eig
      50             :       real :: dx,x0,sig,center,amnr,amni,norm,bf
      51             :       complex :: mmn_arti,amn_arti,amn,psi,wftrial
      52             :       logical :: l_exist,l_miss,l_fe,l_fa,l_fm,l_fm2,l_proj
      53             :       character(len=20) :: fq
      54             :       character(len=32) :: fmt,fmt2
      55             : 
      56           0 :       write(*,*)'create HDWFs input for w90...'
      57             : 
      58             : !      arr_len=3; shifty=0; shiftz=0
      59             : !      if(l_dim(1))then
      60             : !         arr_len=arr_len+1
      61             : !         shifty=shifty+1
      62             : !         shiftz=shiftz+1
      63             : !      endif
      64             : !      if(l_dim(2))then
      65             : !         arr_len=arr_len+1
      66             : !         shiftz=shiftz+1
      67             : !      endif
      68             : !      if(l_dim(3))arr_len=arr_len+1
      69             : 
      70           0 :       call get_dimension(l_dim,arr_len)
      71           0 :       call get_shift(l_dim,shift)
      72             :       
      73           0 :       write(*,*)'dimension:',arr_len
      74           0 :       write(*,*)'x?,y?,z? :',l_dim(1:3)
      75           0 :       if(arr_len.le.3) call juDFT_error("dimension<4",
      76           0 :      >                                  calledby='wann_gwf_commat')
      77             : 
      78           0 :       nkqpts=nkpts*nqpts
      79           0 :       tpi = 2.0*pimach()
      80           0 :       a_arti = latt_const_q
      81           0 :       b_arti = 0.98*latt_const_q
      82             : !      b_arti = 0.5*latt_const_q
      83             : !      b_arti = 0.7*latt_const_q
      84             : !      b_arti = 0.3*latt_const_q
      85             : !      b_arti = 0.05*latt_const_q
      86           0 :       l_fe   = .true.
      87           0 :       l_fa   = .true.
      88           0 :       l_fm   = .true.
      89           0 :       l_fm2  = .true.
      90           0 :       l_exist= .true.
      91             : 
      92             : ! compute auxiliary overlaps and projections
      93           0 :       tau = 1.0*tpi/real(nqpts)/a_arti 
      94           0 :       mmn_arti = 2.0*tpi*tpi*sin(tau*b_arti/2.0)
      95           0 :       mmn_arti = mmn_arti/(tpi*tpi-tau*tau*b_arti*b_arti)
      96           0 :       mmn_arti = mmn_arti/(tau*b_arti)
      97             : c      mmn_arti = 0.5
      98           0 :       write(*,*)'ov(0) =',mmn_arti
      99             : c      CALL wann_gwf_auxovlp(param_vec(:,1),param_vec(:,2),
     100             : c     >                      latt_const_q,mmn_arti)
     101             : 
     102           0 :       amn_arti = cmplx(1.0,0.0)
     103             : 
     104           0 :       open(777,file='ma_aux')
     105           0 :       write(777,'(2(f24.18))')real(mmn_arti),aimag(mmn_arti)
     106           0 :       write(777,'(2(f24.18))')real(amn_arti),aimag(amn_arti)
     107           0 :       write(777,*)a_arti
     108           0 :       write(777,*)b_arti
     109           0 :       close(777)
     110             : 
     111           0 :       bf = b_arti/a_arti
     112           0 :       sig = 1.00!1.2!0.05
     113           0 :       center = 0.0
     114             :       !norm = sqrt(sqrt(2./pimach())/sig)*sqrt(2./bf)
     115           0 :       nx=301*nqpts+1
     116           0 :       dx = 1./real(nx-1)*nqpts
     117           0 :       allocate(xxr(nx),xxi(nx))
     118           0 :       allocate(amn_arti2(nqpts))
     119           0 :       open(555,file='debug_wfaux')
     120           0 :       do q=1,nqpts
     121           0 :        do i=1,nx
     122             :         !x0 = -bf/2. + (i-1)*dx*bf
     123             :         !xx(i) = exp(-x0*x0/sig/sig)*cos(pimach()/bf*x0)
     124           0 :         x0 = -nqpts/2. + (i-1)*dx
     125           0 :         psi=wfpsi(x0,param_vec(3,q),bf)
     126           0 :         wftrial = trialorb(x0,center,sig)
     127           0 :         xxr(i) = real(conjg(wftrial)*psi)
     128           0 :         xxi(i) = aimag(conjg(wftrial)*psi)
     129           0 :         write(555,'(6(f10.6,1x))')param_vec(3,q),x0,
     130           0 :      >       real(psi),aimag(psi),real(wftrial),aimag(wftrial)
     131             :        enddo
     132             :        !call intgz0(xx,dx,nx,amnr,.false.)
     133             :        !amn_arti2(q) = amnr*norm
     134           0 :        call intgz0(xxr,dx,nx,amnr,.false.)
     135           0 :        call intgz0(xxi,dx,nx,amni,.false.)
     136           0 :        amn_arti2(q) = cmplx(amnr,amni)
     137             :       enddo
     138           0 :       close(555)
     139             : 
     140           0 :       open(777,file='debug_amn')
     141           0 :       do q=1,nqpts
     142           0 :        write(777,*)q,amn_arti2(q)
     143             :       enddo
     144           0 :       close(777)
     145           0 :       deallocate(xxr,xxi)
     146             : 
     147             : ! are all necessary files present?
     148           0 :       allocate(feig(nqpts),famn(nqpts),fmmn(nqpts),fmmn2(nqpts))
     149           0 :       l_miss=.false.
     150           0 :       do q=1,nqpts
     151           0 :          write(fq,'(i4.4)')q
     152           0 :          feig(q) = 'WF1_'//trim(fq)//'.eig'    ! energies
     153           0 :          famn(q) = 'WF1_'//trim(fq)//'.amn'    ! projections
     154           0 :          fmmn(q) = 'WF1_'//trim(fq)//'.mmn'    ! k-overlaps
     155           0 :          fmmn2(q)= 'param_'//trim(fq)//'.mmn'  ! q-overlaps
     156             : 
     157           0 :          inquire(file=feig(q) ,exist=l_fe )
     158           0 :          if(l_amn) inquire(file=famn(q) ,exist=l_fa )
     159           0 :          if(l_mmn) inquire(file=fmmn(q) ,exist=l_fm )
     160           0 :          if(l_mmn) inquire(file=fmmn2(q),exist=l_fm2)
     161             : 
     162           0 :          if(.not.l_fe ) write(*,*)'missing: ',feig(q)
     163           0 :          if(.not.l_fa ) write(*,*)'missing: ',famn(q)
     164           0 :          if(.not.l_fm ) write(*,*)'missing: ',fmmn(q)
     165           0 :          if(.not.l_fm2) write(*,*)'missing: ',fmmn2(q)
     166           0 :          if(.not.(l_fe.and.l_fa.and.l_fm.and.l_fm2))l_miss=.true.
     167             :       enddo
     168           0 :       if(l_mmn) inquire(file='bkqpts',exist=l_exist)
     169           0 :       if(.not.l_exist) write(*,*)'missing: bkqpts'
     170           0 :       inquire(file='proj',exist=l_proj)
     171           0 :       if(.not.l_proj) write(*,*)'missing: proj'
     172           0 :       if(l_miss.or.(.not.l_exist).or.(.not.l_proj))
     173           0 :      >   call juDFT_error("missing file(s) for HDWFs")
     174             : 
     175             : ! get number of bands and wfs from proj
     176           0 :       open(405,file='proj',status='old')
     177           0 :       read(405,*)nwfs,nbnd
     178           0 :       close(405)
     179           0 :       write(*,*)'nbnd=',nbnd
     180           0 :       write(*,*)'nwfs=',nwfs
     181             : 
     182             : 
     183             : c*****************************c
     184             : c     COMPOSITE .EIG FILE     c
     185             : c*****************************c
     186           0 :       open(305,file='WF1_gwf.eig')
     187           0 :       do q=1,nqpts
     188           0 :          open(405,file=feig(q))
     189           0 :          do k=1,nkpts
     190           0 :             ikqpt=get_index_kq(k,q,nkpts)
     191           0 :             do i=1,nbnd
     192           0 :                read(405,*)b,nk,eig
     193           0 :                write(305,'(2i12,f19.13)')b,ikqpt,eig
     194             :             enddo
     195             :          enddo
     196           0 :          close(405)!,status='delete')
     197             :       enddo
     198           0 :       close(305)
     199             : 
     200             : 
     201             : c*****************************c
     202             : c     COMPOSITE .AMN FILE     c
     203             : c*****************************c
     204           0 :       if(.not.l_amn) goto 100 ! skip amn part
     205             : 
     206           0 :       if(l_unformatted) then
     207           0 :          allocate(tmp_a(nbnd,nwfs,nkqpts))
     208           0 :          do q=1,nqpts
     209           0 :             is = get_index_kq(1,q,nkpts)
     210           0 :             ie = get_index_kq(nkpts,q,nkpts)
     211           0 :             open(405,file=famn(q),form='unformatted')
     212           0 :             read(405)b,k,nwf
     213           0 :             read(405)tmp_a(:,:,is:ie)
     214           0 :             do i=1,nwfs
     215           0 :             tmp_a(:,i,is:ie) = tmp_a(:,i,is:ie)*amn_arti2(q)
     216             :             enddo
     217           0 :             close(405)!,status='delete')
     218             :          enddo
     219             : 
     220           0 :          open(306,file='WF1_gwf.amn',form='unformatted')
     221           0 :          write(306)nbnd,nkqpts,nwfs
     222           0 :          write(306)tmp_a!*amn_arti
     223           0 :          close(306)
     224           0 :          deallocate(tmp_a)
     225             :       else
     226           0 :          open(306,file='WF1_gwf.amn')
     227           0 :          write(306,*)'Projections for HDWFs'
     228           0 :          write(306,'(i5,i7,i5)')nbnd,nkqpts,nwfs
     229             :          
     230           0 :          do q=1,nqpts
     231           0 :            open(405,file=famn(q))
     232           0 :            read(405,*)!title
     233           0 :            read(405,*)b,k,nwf
     234           0 :            do k=1,nkpts
     235           0 :             ikqpt=get_index_kq(k,q,nkpts)
     236           0 :             do nwf=1,nwfs
     237           0 :              do i=1,nbnd
     238           0 :               read(405,*)b,n,nk,mmn_r,mmn_i
     239           0 :               amn = cmplx(mmn_r,mmn_i)
     240           0 :               amn = amn*amn_arti2(q) !*amn_arti
     241             :               write(306,'(i5,i5,i7,3x,2f18.12)')
     242           0 :      >                     b,n,ikqpt,real(amn),aimag(amn)
     243             :              enddo
     244             :             enddo
     245             :            enddo
     246           0 :            close(405,status='delete')
     247             :          enddo
     248           0 :          close(306)
     249             :       endif!l_unformatted
     250             : 
     251             :  100  continue
     252             : 
     253             : 
     254             : c*****************************c
     255             : c     COMPOSITE .MMN FILE     c
     256             : c*****************************c
     257           0 :       if(.not.l_mmn) goto 200 ! skip mmn part
     258             : 
     259           0 :       write(fmt,'(a,i1,a)')'(2i6,3x,',arr_len,'i4)'
     260           0 :       write(fmt2,'(a,i1,a)')'(i7,i7,3x,',arr_len,'i4)'
     261             : 
     262           0 :       open (202,file='bkqpts',form='formatted',status='old')
     263           0 :       rewind (202)
     264           0 :       read (202,'(i4)') nntot_kq
     265           0 :       write (*,*) 'nntot_kq=',nntot_kq
     266           0 :       allocate ( gb_kq(arr_len,nntot_kq,nkqpts),
     267           0 :      &           bpt_kq(nntot_kq,nkqpts))
     268           0 :       do ikqpt=1,nkqpts
     269           0 :        do nn=1,nntot_kq
     270             :           read (202,fmt)
     271           0 :      &      ikqpt_help,bpt_kq(nn,ikqpt),(gb_kq(i,nn,ikqpt),i=1,arr_len)
     272             :        enddo
     273             :       enddo
     274           0 :       close (202)
     275             : 
     276           0 :       if(l_unformatted) then
     277           0 :          allocate(mmnk(nbnd,nbnd,nntot,nkpts))
     278           0 :          allocate(mmnq(nbnd,nbnd,nntot_q,nkpts))
     279           0 :          allocate(g1(3,nntot,nkpts),g2(3,nntot_q))
     280           0 :          allocate(b1(nntot,nkpts),b2(nntot_q))
     281           0 :          allocate(mmn(nbnd,nbnd,nntot_kq,nkqpts))
     282             : 
     283           0 :          do q=1,nqpts
     284           0 :            is = get_index_kq(1,q,nkpts)
     285           0 :            ie = get_index_kq(nkpts,q,nkpts)
     286             :          
     287           0 :            open(405,file=fmmn(q),form='unformatted')
     288           0 :            read(405)b,k,nn
     289           0 :            read(405)b1,g1
     290           0 :            read(405)mmnk
     291           0 :            close(405,status='delete')
     292             : 
     293           0 :            open(405,file=fmmn2(q),form='unformatted')
     294           0 :            read(405)b,k,nn
     295           0 :            read(405)b2,g2
     296           0 :            read(405)mmnq
     297           0 :            close(405,status='delete')
     298           0 :            mmnq = mmnq*conjg(mmn_arti)
     299             : 
     300           0 :            do k=1,nkpts
     301           0 :              ikqpt=get_index_kq(k,q,nkpts)           
     302           0 :              do nn=1,nntot_kq
     303           0 :                kqb = bpt_kq(nn,ikqpt)
     304           0 :                kb=get_index_k(kqb,nkpts)
     305           0 :                qb=get_index_q(kqb,nkpts)
     306           0 :                if(q.eq.qb) then ! k-overlap
     307             :                   nk=get_index_nn_k(bpt(:,k),nntot,kb,
     308             :      >                         gb_kq(:,nn,ikqpt),
     309           0 :      >                         gb(1:3,1:nntot,k),arr_len)
     310           0 :                   mmn(:,:,nn,ikqpt)=mmnk(:,:,nk,k)              
     311             :                   !write(*,*)'k neig',k,kb,nk
     312           0 :                elseif(k.eq.kb) then ! q-overlap
     313             :                   nq=get_index_nn_q(bpt_q(:,q),nntot_q,qb,
     314             :      >                         gb_kq(:,nn,ikqpt),
     315             :      >                         gb_q(1:3,1:nntot_q,q),
     316           0 :      >                         arr_len,shift,l_dim)
     317           0 :                   mmn(:,:,nn,ikqpt)=mmnq(:,:,nq,k)
     318             :                   !write(*,*)'q neig',q,qb,nq
     319             :                else ! otherwise problem
     320             :                   call juDFT_error("problem in overlap mmn gwf",
     321           0 :      >                             calledby="wann_gwf_commat")
     322             :                endif
     323             :              enddo!nn
     324             :            enddo!k
     325             :          enddo!q
     326           0 :          deallocate(b1,b2,g1,g2,mmnk,mmnq)
     327             : 
     328           0 :          open(307,file='WF1_gwf.mmn',form='unformatted')
     329           0 :          write(307)nbnd,nkqpts,nntot_kq
     330           0 :          write(307)bpt_kq,gb_kq
     331           0 :          write(307)mmn
     332           0 :          close(307)
     333           0 :          deallocate(mmn)
     334             :       else
     335           0 :          allocate(mk(nbnd,nbnd,nntot))
     336           0 :          allocate(mq(nbnd,nbnd,nntot_q))
     337           0 :          open(307,file='WF1_gwf.mmn')
     338           0 :          write(307,*)'Overlaps for HDWFs'
     339           0 :          write(307,'(i5,i7,i5)')nbnd,nkqpts,nntot_kq
     340             : 
     341           0 :          do q=1,nqpts       
     342           0 :            open(405,file=fmmn(q))
     343           0 :            read(405,*)!title
     344           0 :            read(405,*)b,k,nn
     345           0 :            open(406,file=fmmn2(q))
     346           0 :            read(406,*)!title
     347           0 :            read(406,*)b,k,nn
     348           0 :            do k=1,nkpts
     349           0 :              ikqpt=get_index_kq(k,q,nkpts)           
     350             : 
     351           0 :              do nn=1,nntot ! read k-overlaps
     352           0 :                kb = bpt(nn,k)
     353           0 :                read(405,*)nkp,b,g(1:3)
     354           0 :                do i=1,nbnd
     355           0 :                   do j=1,nbnd
     356           0 :                      read(405,*)mmn_r,mmn_i
     357           0 :                      mk(j,i,nn)=cmplx(mmn_r,mmn_i)
     358             :                   enddo
     359             :                enddo               
     360             :              enddo  
     361             : 
     362           0 :              do nn=1,nntot_q ! read q-overlaps
     363           0 :                qb = bpt_q(nn,q)
     364           0 :                read(406,*)nkp,b,g(1:3)
     365           0 :                do i=1,nbnd
     366           0 :                   do j=1,nbnd
     367           0 :                      read(406,*)mmn_r,mmn_i
     368           0 :                      mq(j,i,nn)=cmplx(mmn_r,mmn_i)
     369             :                   enddo
     370             :                enddo                              
     371             :              enddo
     372           0 :              mq = mq*conjg(mmn_arti)   
     373             : 
     374           0 :              do nn=1,nntot_kq ! write composite overlaps
     375           0 :                kqb = bpt_kq(nn,ikqpt)
     376           0 :                kb=get_index_k(kqb,nkpts)
     377           0 :                qb=get_index_q(kqb,nkpts)
     378           0 :                write(307,fmt2)ikqpt,kqb,gb_kq(1:arr_len,nn,ikqpt)
     379           0 :                if(q.eq.qb) then ! k-overlap
     380             :                   nk=get_index_nn_k(bpt(:,k),nntot,kb,
     381             :      >                         gb_kq(:,nn,ikqpt),
     382           0 :      >                         gb(1:3,1:nntot,k),arr_len)
     383           0 :                   do i=1,nbnd
     384           0 :                      do j=1,nbnd
     385             :                         write(307,'(2f24.18)')
     386           0 :      >                      real(mk(j,i,nk)),aimag(mk(j,i,nk))
     387             :                      enddo
     388             :                   enddo
     389           0 :                elseif(k.eq.kb) then ! q-overlap
     390             :                   nq=get_index_nn_q(bpt_q(:,q),nntot_q,qb,
     391             :      >                         gb_kq(:,nn,ikqpt),
     392             :      >                         gb_q(1:3,1:nntot_q,q),
     393           0 :      >                         arr_len,shift,l_dim)
     394           0 :                   do i=1,nbnd
     395           0 :                      do j=1,nbnd
     396             :                         write(307,'(2f24.18)')
     397           0 :      >                      real(mq(j,i,nq)),aimag(mq(j,i,nq))
     398             :                      enddo
     399             :                   enddo
     400             :                else ! otherwise problem
     401             :                   call juDFT_error("problem in overlap mmn gwf",
     402           0 :      >                             calledby="wann_gwf_commat")
     403             :                endif
     404             :              enddo!nn
     405             :            enddo!k
     406           0 :            close(406,status='delete')
     407           0 :            close(405,status='delete')
     408             :          enddo!q           
     409           0 :          close(307)
     410           0 :          deallocate(mk,mq)
     411             :       endif!l_unformatted
     412             : 
     413           0 :       deallocate(gb_kq,bpt_kq)
     414             :  200  continue
     415           0 :       deallocate(feig,famn,fmmn,fmmn2)
     416           0 :       deallocate(amn_arti2)
     417             : 
     418             :       contains
     419             : 
     420             : 
     421           0 :       complex function wfu(x0,qpt,bb)
     422             :       implicit none
     423             :       real,intent(in) :: x0,qpt,bb
     424             :       integer :: i0
     425             :       real :: qx,xshift
     426             :       complex :: pha
     427             : 
     428           0 :       wfu = cmplx(0.,0.)
     429           0 :       i0=floor(x0+0.5)
     430           0 :       xshift=x0-real(i0)
     431           0 :       if(abs(xshift) .ge. bb/2.) return
     432             : 
     433           0 :       qx = 2.*pimach()*qpt*xshift
     434           0 :       pha = cmplx( cos(qx), -sin(qx) )
     435           0 :       wfu = pha*sqrt(2./bb)*cos(pimach()*xshift/bb)
     436           0 :       end function wfu
     437             : 
     438           0 :       complex function wfpsi(x0,qpt,bb)
     439             :       implicit none
     440             :       real,intent(in) :: x0,qpt,bb
     441             :       complex :: pha
     442             :       real :: qx
     443             :  
     444           0 :       qx = 2.*pimach()*qpt*x0
     445           0 :       pha = cmplx( cos(qx), sin(qx) )
     446           0 :       wfpsi = pha*wfu(x0,qpt,bb)      
     447           0 :       end function wfpsi
     448             : 
     449             :       
     450           0 :       real function trialorb(x0,center,sigma)
     451             :       implicit none
     452             :       real,intent(in) :: x0,center,sigma
     453             :       real :: xp
     454             : 
     455           0 :       xp = (x0-center)/sigma
     456           0 :       trialorb = exp(-xp*xp)*sqrt(sqrt(2./pimach())/sigma)
     457           0 :       end function trialorb
     458             : 
     459             :       end subroutine wann_gwf_commat    
     460             :       end module m_wann_gwf_commat

Generated by: LCOV version 1.13