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

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : 
       7             :       module m_wann_projmethod
       8             :       use m_juDFT
       9             :       contains
      10           0 :       subroutine wann_projmethod(
      11             :      >               fullnkpts,
      12             :      >               l_projmethod,l_bestproj,
      13           0 :      >               ikpt,nwfs,nslibd,amn,eig,
      14           0 :      <               psiw,hwfr)
      15             : c**********************************************************************
      16             : c     Construct Wannier functions from the projections of the
      17             : c     Bloch functions onto the trial orbitals. These projections
      18             : c     are provided as input to this subroutine in the amn-matrix.
      19             : c     The Wannier functions are orthonormalized. The orthonormalization
      20             : c     requires the construction of square root of the overlap 
      21             : c     matrix (omn). Two different algorithms for the calculation of the
      22             : c     square root are available (projmethod/bestproj).
      23             : c
      24             : c     YM && FF
      25             : c**********************************************************************
      26             :       use m_types
      27             :       implicit none
      28             :       integer, intent(in)    :: fullnkpts
      29             :       logical, intent(in)    :: l_projmethod
      30             :       logical, intent(in)    :: l_bestproj
      31             :       integer, intent(in)    :: ikpt
      32             :       integer, intent(in)    :: nwfs
      33             :       integer, intent(in)    :: nslibd
      34             :       complex, intent(in)    :: amn(:,:,:)
      35             :       real,    intent(in)    :: eig(:)
      36             : 
      37             :       complex, intent(inout) :: psiw(:,:,:)
      38             :       complex, intent(inout) :: hwfr(:,:)
      39             : 
      40           0 :       real                   :: ei(nwfs)
      41           0 :       complex, allocatable   :: work(:)
      42           0 :       integer, allocatable   :: iwork(:)
      43           0 :       real, allocatable      :: rwork(:)
      44             :       integer                :: info,lwork,lrwork,liwork,lee
      45             :       integer                :: nwf,nwfp,i,n,lda,j
      46           0 :       complex                :: omn(nwfs,nwfs),vec(nwfs,nwfs)
      47           0 :       complex                :: smn(nwfs,nwfs)
      48           0 :       complex, allocatable   :: amn_copy(:,:)
      49             :       character              :: jobu*1,jobv*1
      50             :       character              :: jobz*1,uplo*1
      51           0 :       complex, allocatable   :: sunit(:,:)
      52           0 :       real,    allocatable   :: sdiag(:)
      53           0 :       complex, allocatable   :: umn(:,:),vmn(:,:)
      54             :       integer                :: ldu,ldv,m
      55           0 :       complex                :: hwfk(nwfs,nwfs)
      56             : 
      57           0 :       if(l_projmethod) then
      58           0 :          omn=cmplx(0.0,0.0)
      59             : c..      now we determine the Omn matrix, PRB 71,125119, Eq.(13)
      60           0 :          do nwf = 1,nwfs
      61           0 :           do nwfp = 1,nwfs
      62           0 :            do i = 1,nslibd
      63             :               omn(nwf,nwfp) = omn(nwf,nwfp) +
      64           0 :      +               conjg(amn(i,nwf,ikpt))*amn(i,nwfp,ikpt)
      65             :            enddo
      66             :           enddo
      67             :          enddo
      68             : 
      69             : c        lwork = 2*nwfs + nwfs*nwfs
      70             : c        lrwork = 1 + 4*nwfs + 5*nwfs*nwfs
      71             : c        liwork = 2 + 6*nwfs
      72           0 :          lee = log( dble(nwfs) )/log(2.d0) + 1
      73           0 :          lwork = 1 + 5*nwfs + 2*nwfs*lee + 3*(nwfs**2)
      74           0 :          lrwork = 1 + 4*nwfs + 2*nwfs*lee + 3*(nwfs**2)
      75           0 :          liwork = 2 + 5*nwfs +1
      76             : 
      77           0 :          allocate (work(lwork),rwork(lrwork),iwork(liwork))
      78             : 
      79           0 :          jobz = 'V' ; uplo = 'L' ; n = nwfs ; lda = nwfs
      80             : 
      81             :          call zheevd(jobz,uplo,n,vec,lda,ei,work,lwork,
      82           0 :      &             rwork,lrwork,iwork,liwork,info)
      83             : 
      84           0 :          if (info.lt.0) write (6,*)
      85           0 :      &               'ith argument had an illegal value ',info
      86           0 :          IF (info>0)  CALL juDFT_error("not converged diagonalization"
      87           0 :      +        ,calledby ="wann_projmethod")
      88             : 
      89           0 :          do i = 1,nwfs
      90           0 :           if (ei(i).le.(-1.e-6)) then
      91           0 :            print*,'eigenvalue is less than zero'
      92           0 :           elseif (abs(ei(i)).le.1.e-16) then
      93           0 :            print*,'eigenvalue is very close to zero:',ei(i)
      94             :           endif
      95             :          enddo
      96             : 
      97             : c..      constructing the square root of the O matrix, Eq.(14)
      98             : 
      99           0 :          smn(:,:) = cmplx(0.,0.)
     100             : 
     101           0 :          do i = 1,nwfs
     102           0 :           do j = 1,nwfs
     103           0 :            do n = 1,nwfs
     104             :             smn(i,j) = smn(i,j) +
     105           0 :      +         conjg(vec(i,n))*(1./sqrt(ei(n)))*vec(j,n)
     106             :            enddo
     107             :           enddo
     108             :          enddo
     109             : 
     110           0 :          deallocate (work,rwork,iwork)
     111             : c..      now constructing the overlaps of the wavefunctions
     112             : c..      with the wannier functions in reciprocal space, Eq.(16)
     113             : 
     114           0 :          psiw(:,:,ikpt) = cmplx(0.,0.)
     115           0 :          do i = 1,nslibd
     116           0 :           do j = 1,nwfs
     117           0 :            do n = 1,nwfs
     118             :                psiw(i,j,ikpt) = psiw(i,j,ikpt) +
     119           0 :      +                       smn(j,n)*amn(i,n,ikpt)
     120             :            enddo
     121             :           enddo
     122             :          enddo
     123             :       endif !wann%l_projmethod
     124             : 
     125           0 :       if(l_bestproj) then
     126             : 
     127             : c..      psiw matrix has the meaning of the Umn rotational matrix
     128             : c..      applied to the wavefunctions at each point in the BZone
     129             : c..      here it is calculated differently,
     130             : c..      according to the PRB 65,035109 (Eq.23)
     131             : c..
     132             : c..      the singe value decomposition of the Amn matrix is found:
     133             : c..                      A = U*SS*V
     134             : c..      Then the psiw matrix, psiw = Amn*Smn is given by:
     135             : c..                     psiw = U*1*V
     136             : 
     137           0 :          allocate (amn_copy(nslibd,nwfs),sunit(nslibd,nwfs))
     138           0 :          allocate (umn(nslibd,nslibd),vmn(nwfs,nwfs))
     139           0 :          allocate (sdiag(max(nslibd,nwfs)))
     140             : 
     141             :          lwork = max( 3*min(nslibd,nwfs) + max(nslibd,nwfs),
     142           0 :      &                        5*min(nslibd,nwfs) )
     143             : 
     144           0 :          allocate ( work(lwork),rwork(max(1,5*min(nslibd,nwfs))) )
     145             : !
     146             : !        Compute the eigenvalues and eigenvectors.
     147             : !
     148           0 :          jobu = 'A' ; jobv = 'A'
     149           0 :          lda = nslibd ; ldu = nslibd ; ldv = nwfs
     150             : !
     151             : !        The input matrix is destroyed by the routine.  Since we need to keep
     152             : !        it around, we only pass a copy to the routine.
     153             : !
     154           0 :          amn_copy(1:nslibd,1:nwfs) = amn(1:nslibd,1:nwfs,ikpt)
     155             : 
     156             :          call zgesvd(jobu,jobv,nslibd,nwfs,amn_copy,lda,
     157           0 :      >             sdiag,umn,ldu,vmn,ldv,work,lwork,rwork,info)
     158             : 
     159           0 :          if (info /= 0) then
     160           0 :             write (*,*) 'LAPACK routine ZGESVD returned a nonzero'
     161           0 :             write (*,*) 'value of the error flag, INFO =',info
     162             :          end if
     163             : !
     164             : !        Make the MxN identical Sunit matrix
     165             : !
     166           0 :          sunit(1:nslibd,1:nwfs) = cmplx(0.,0.)
     167           0 :          do i = 1,min(nslibd,nwfs)
     168           0 :             sunit(i,i) = cmplx(1.,0.)
     169             :          end do
     170             : 
     171             : 
     172             : !        and finally the psiw matrix
     173             : 
     174           0 :          psiw(:,:,ikpt) = cmplx(0.,0.)
     175           0 :          do i = 1,nslibd
     176           0 :           do j = 1,nwfs
     177           0 :            do n = 1,nslibd
     178           0 :             do m = 1,nwfs
     179             :                psiw(i,j,ikpt) = psiw(i,j,ikpt) +
     180           0 :      +                     umn(i,n)*sunit(n,m)*vmn(m,j)
     181             :             enddo
     182             :            enddo
     183             :           enddo
     184             :          enddo
     185             : 
     186           0 :          deallocate (work,rwork,amn_copy,sunit,vmn,umn,sdiag)
     187             : 
     188           0 :          write (*,*) 'The Psiw matrix was calculated via SVD'
     189             : 
     190             :       endif   !wann%l_bestproj
     191             : 
     192             : 
     193             : 
     194             : c..constructing the WF-hamiltonian in reciprocal space Eq.(23)
     195             : 
     196           0 :       hwfk(:,:) = cmplx(0.,0.)
     197             : 
     198           0 :       do i = 1,nwfs
     199           0 :          do j = 1,nwfs
     200           0 :             do n = 1,nslibd
     201             :                hwfk(i,j) = hwfk(i,j) +
     202           0 :      +            eig(n)*psiw(n,i,ikpt)*conjg(psiw(n,j,ikpt))
     203             : 
     204             :             enddo
     205             :          enddo
     206             :       enddo
     207             : 
     208             : c..   adding up the k-point reciprocal hamiltonians Eq.(20)
     209             : c..   to get the hopping elements inside the same unit cell
     210             : 
     211           0 :       hwfr=hwfr+hwfk/fullnkpts
     212             : 
     213             : 
     214             : c..   now we diagonalize the reciprocal hamiltonian for the
     215             : c..   purpose of testing
     216             : 
     217           0 :       do nwf = 1,nwfs
     218           0 :        do nwfp = 1,nwfs
     219           0 :         vec(nwf,nwfp) = hwfk(nwf,nwfp)
     220             :        enddo
     221             :       enddo
     222             : 
     223           0 :       lee = log( dble(nwfs) )/log(2.d0) + 1
     224           0 :       lwork = 1 + 5*nwfs + 2*nwfs*lee + 3*(nwfs**2)
     225           0 :       lrwork = 1 + 4*nwfs + 2*nwfs*lee + 3*(nwfs**2)
     226           0 :       liwork = 2 + 5*nwfs +1
     227             : 
     228           0 :       allocate (work(lwork),rwork(lrwork),iwork(liwork))
     229             : 
     230           0 :       jobz = 'V' ; uplo = 'L' ; n = nwfs ; lda = nwfs
     231             : 
     232             :       call zheevd(jobz,uplo,n,vec,lda,ei,work,lwork,
     233           0 :      &             rwork,lrwork,iwork,liwork,info)
     234             : 
     235           0 :       if (info.lt.0) write (6,*)
     236           0 :      &               'ith argument had an illegal value ',info
     237           0 :       IF (info>0)  CALL juDFT_error("not converged diagonalization"
     238           0 :      +     ,calledby ="wann_projmethod")
     239             : 
     240           0 :       do i = 1,nwfs
     241           0 :          write(6,*) ei(i)
     242             :       enddo
     243             : 
     244           0 :       deallocate (work,rwork,iwork)
     245             : 
     246             : 
     247           0 :       end subroutine wann_projmethod
     248             :       end module m_wann_projmethod

Generated by: LCOV version 1.13