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

Generated by: LCOV version 2.0-1