src/seceqnhf.f90

   1
   2 ! Copyright (C) 2006 J. K. Dewhurst and S. Sharma.
   3 ! This file is distributed under the terms of the GNU General Public License.
   4 ! See the file COPYING for license details.
   5
   6 subroutine seceqnhf(ikp,evecsvp)
   7 use modmain
   8 implicit none
   9 ! arguments
  10 integer, intent(in) :: ikp
  11 complex(8), intent(inout) :: evecsvp(nstsv,nstsv)
  12 ! local variables
  13 integer is,ia,ias,ir,irc
  14 integer ngknr,ik,jk,ist1,ist2,ist3
  15 integer iq,ig,iv(3),igq0
  16 integer lmax,lwork,info
  17 real(8) cfq,v(3),t1
  18 complex(8) zrho01,zrho02,zt1,zt2
  19 ! allocatable arrays
  20 integer, allocatable :: igkignr(:)
  21 real(8), allocatable :: vgklnr(:,:)
  22 real(8), allocatable :: vgkcnr(:,:)
  23 real(8), allocatable :: gkcnr(:)
  24 real(8), allocatable :: tpgkcnr(:,:)
  25 real(8), allocatable :: vgqc(:,:)
  26 real(8), allocatable :: tpgqc(:,:)
  27 real(8), allocatable :: gqc(:)
  28 real(8), allocatable :: jlgqr(:,:,:)
  29 real(8), allocatable :: evalsvp(:)
  30 real(8), allocatable :: evalsvnr(:)
  31 real(8), allocatable :: rwork(:)
  32 real(8), allocatable :: rfmt(:,:,:)
  33 complex(8), allocatable :: h(:,:)
  34 complex(8), allocatable :: vmat(:,:)
  35 complex(8), allocatable :: apwalm(:,:,:,:)
  36 complex(8), allocatable :: evecfv(:,:)
  37 complex(8), allocatable :: evecsv(:,:)
  38 complex(8), allocatable :: sfacgknr(:,:)
  39 complex(8), allocatable :: ylmgq(:,:)
  40 complex(8), allocatable :: sfacgq(:,:)
  41 complex(8), allocatable :: wfmt1(:,:,:,:,:)
  42 complex(8), allocatable :: wfmt2(:,:,:,:,:)
  43 complex(8), allocatable :: wfir1(:,:,:)
  44 complex(8), allocatable :: wfir2(:,:,:)
  45 complex(8), allocatable :: zrhomt(:,:,:)
  46 complex(8), allocatable :: zrhoir(:)
  47 complex(8), allocatable :: zpchg(:)
  48 complex(8), allocatable :: zvclmt(:,:,:)
  49 complex(8), allocatable :: zvclir(:)
  50 complex(8), allocatable :: zfmt(:,:)
  51 complex(8), allocatable :: work(:)
  52 ! external functions
  53 complex(8) zfinp,zfmtinp
  54 external zfinp,zfmtinp
  55 !$OMP CRITICAL
  56 write(*,'("Info(seceqnhf): ",I6," of ",I6," k-points")') ikp,nkpt
  57 !$OMP END CRITICAL
  58 ! allocate local arrays
  59 allocate(igkignr(ngkmax))
  60 allocate(vgklnr(3,ngkmax))
  61 allocate(vgkcnr(3,ngkmax))
  62 allocate(gkcnr(ngkmax))
  63 allocate(tpgkcnr(2,ngkmax))
  64 allocate(vgqc(3,ngvec))
  65 allocate(tpgqc(2,ngvec))
  66 allocate(gqc(ngvec))
  67 allocate(jlgqr(0:lmaxvr+npsden+1,ngvec,nspecies))
  68 allocate(evalsvp(nstsv))
  69 allocate(evalsvnr(nstsv))
  70 allocate(rwork(3*nstsv))
  71 allocate(h(nstsv,nstsv))
  72 allocate(vmat(nstsv,nstsv))
  73 allocate(apwalm(ngkmax,apwordmax,lmmaxapw,natmtot))
  74 allocate(evecfv(nmatmax,nstfv))
  75 allocate(evecsv(nstsv,nstsv))
  76 allocate(sfacgknr(ngkmax,natmtot))
  77 allocate(ylmgq(lmmaxvr,ngvec))
  78 allocate(sfacgq(ngvec,natmtot))
  79 allocate(wfmt1(lmmaxvr,nrcmtmax,natmtot,nspinor,nstsv))
  80 allocate(wfmt2(lmmaxvr,nrcmtmax,natmtot,nspinor,nstsv))
  81 allocate(wfir1(ngrtot,nspinor,nstsv))
  82 allocate(wfir2(ngrtot,nspinor,nstsv))
  83 allocate(zrhomt(lmmaxvr,nrcmtmax,natmtot))
  84 allocate(zrhoir(ngrtot))
  85 allocate(zpchg(natmtot))
  86 allocate(zvclmt(lmmaxvr,nrcmtmax,natmtot))
  87 allocate(zvclir(ngrtot))
  88 allocate(zfmt(lmmaxvr,nrcmtmax))
  89 lwork=2*nstsv
  90 allocate(work(lwork))
  91 allocate(rfmt(lmmaxvr,nrcmtmax,natmtot))
  92 ! coefficient of long-range term
  93 cfq=0.5d0*(omega/pi)**2
  94 ! set the point charges to zero
  95 zpchg(:)=0.d0
  96 ! get the eigenvalues/vectors from file for input k-point
  97 call getevalsv(vkl(1,ikp),evalsvp)
  98 call getevecfv(vkl(1,ikp),vgkl(1,1,ikp,1),evecfv)
  99 ! find the matching coefficients
 100 call match(ngk(ikp,1),gkc(1,ikp,1),tpgkc(1,1,ikp,1),sfacgk(1,1,ikp,1),apwalm)
 101 ! calculate the wavefunctions for all states for the input k-point
 102 call genwfsv(.false.,ngk(ikp,1),igkig(1,ikp,1),evalsvp,apwalm,evecfv,evecsvp, &
 103  wfmt1,wfir1)
 104 ! compute the new kinetic matrix elements
 105 call zgemm('N','N',nstsv,nstsv,nstsv,zone,kinmatc(1,1,ikp),nstsv,evecsvp, &
 106  nstsv,zzero,vmat,nstsv)
 107 call zgemm('C','N',nstsv,nstsv,nstsv,zone,evecsvp,nstsv,vmat,nstsv,zzero,h, &
 108  nstsv)
 109 ! convert muffin-tin Coulomb potential to spherical coordinates
 110 do is=1,nspecies
 111   do ia=1,natoms(is)
 112     ias=idxas(ia,is)
 113     irc=0
 114     do ir=1,nrmt(is),lradstp
 115       irc=irc+1
 116       call dgemv('N',lmmaxvr,lmmaxvr,1.d0,rbshtapw,lmmaxapw,vclmt(1,ir,ias), &
 117        1,0.d0,rfmt(1,irc,ias),1)
 118     end do
 119   end do
 120 end do
 121 ! compute the Coulomb matrix elements and add
 122 call genvmatk(rfmt,vclir,wfmt1,wfir1,vmat)
 123 h(:,:)=h(:,:)+vmat(:,:)
 124 ! zero the non-local matrix elements for passed k-point
 125 vmat(:,:)=0.d0
 126 ! start loop over non-reduced k-point set
 127 do ik=1,nkptnr
 128 ! find the equivalent reduced k-point
 129   iv(:)=ivknr(:,ik)
 130   jk=ikmap(iv(1),iv(2),iv(3))
 131 ! generate G+k vectors
 132   call gengpvec(vklnr(1,ik),vkcnr(1,ik),ngknr,igkignr,vgklnr,vgkcnr,gkcnr, &
 133    tpgkcnr)
 134 ! get the eigenvalues/vectors from file for non-reduced k-point
 135   call getevalsv(vklnr(1,ik),evalsvnr)
 136   call getevecfv(vklnr(1,ik),vgklnr,evecfv)
 137   call getevecsv(vklnr(1,ik),evecsv)
 138 ! generate the structure factors
 139   call gensfacgp(ngknr,vgkcnr,ngkmax,sfacgknr)
 140 ! find the matching coefficients
 141   call match(ngknr,gkcnr,tpgkcnr,sfacgknr,apwalm)
 142 ! determine q-vector
 143   iv(:)=ivk(:,ikp)-ivknr(:,ik)
 144   iv(:)=modulo(iv(:),ngridk(:))
 145   iq=iqmap(iv(1),iv(2),iv(3))
 146   v(:)=vkc(:,ikp)-vkcnr(:,ik)
 147   do ig=1,ngvec
 148 ! determine G+q vectors
 149     vgqc(:,ig)=vgc(:,ig)+v(:)
 150 ! G+q-vector length and (theta, phi) coordinates
 151     call sphcrd(vgqc(1,ig),gqc(ig),tpgqc(1,ig))
 152 ! spherical harmonics for G+q-vector
 153     call genylm(lmaxvr,tpgqc(1,ig),ylmgq(1,ig))
 154   end do
 155 ! structure factors for G+q
 156   call gensfacgp(ngvec,vgqc,ngvec,sfacgq)
 157 ! find the shortest G+q-vector
 158   call findigp0(ngvec,gqc,igq0)
 159 ! compute the required spherical Bessel functions
 160   lmax=lmaxvr+npsden+1
 161   call genjlgpr(lmax,gqc,jlgqr)
 162 ! calculate the wavefunctions for all states
 163   call genwfsv(.false.,ngknr,igkignr,evalsvnr,apwalm,evecfv,evecsv,wfmt2,wfir2)
 164   do ist3=1,nstsv
 165     if (occsv(ist3,jk).gt.epsocc) then
 166       do ist2=1,nstsv
 167 ! calculate the complex overlap density
 168         call vnlrho(.true.,wfmt2(1,1,1,1,ist3),wfmt1(1,1,1,1,ist2), &
 169          wfir2(1,1,ist3),wfir1(1,1,ist2),zrhomt,zrhoir)
 170 ! calculate the Coulomb potential
 171         call zpotcoul(nrcmt,nrcmtmax,nrcmtmax,rcmt,igq0,gqc,jlgqr,ylmgq, &
 172          sfacgq,zpchg,zrhomt,zrhoir,zvclmt,zvclir,zrho02)
 173 !----------------------------------------------!
 174 !     valence-valence-valence contribution     !
 175 !----------------------------------------------!
 176         do ist1=1,ist2
 177 ! calculate the complex overlap density
 178           call vnlrho(.true.,wfmt2(1,1,1,1,ist3),wfmt1(1,1,1,1,ist1), &
 179            wfir2(1,1,ist3),wfir1(1,1,ist1),zrhomt,zrhoir)
 180           zt1=zfinp(.true.,zrhomt,zvclmt,zrhoir,zvclir)
 181 ! compute the density coefficient of the smallest G+q-vector
 182           call zrhoqint(gqc(igq0),ylmgq(1,igq0),ngvec,sfacgq(igq0,1),zrhomt, &
 183            zrhoir,zrho01)
 184           zt2=cfq*wiq2(iq)*(conjg(zrho01)*zrho02)
 185           t1=occsv(ist3,jk)/occmax
 186           vmat(ist1,ist2)=vmat(ist1,ist2)-t1*(wkptnr(ik)*zt1+zt2)
 187         end do
 188       end do
 189     end if
 190   end do
 191 ! end loop over non-reduced k-point set
 192 end do
 193 ! add the non-local matrix elements to Hamiltonian
 194 h(:,:)=h(:,:)+vmat(:,:)
 195 ! diagonalise the Hartree-Fock Hamiltonian (eigenvalues in global array)
 196 call zheev('V','U',nstsv,h,nstsv,evalsv(1,ikp),work,lwork,rwork,info)
 197 if (info.ne.0) then
 198   write(*,*)
 199   write(*,'("Error(seceqnhf): diagonalisation of the Hartree-Fock Hamiltonian &
 200    &failed")')
 201   write(*,'(" for k-point ",I8)') ikp
 202   write(*,'(" ZHEEV returned INFO = ",I8)') info
 203   write(*,*)
 204   stop
 205 end if
 206 ! apply unitary transformation to second-variational states
 207 evecsv(:,:)=evecsvp(:,:)
 208 call zgemm('N','N',nstsv,nstsv,nstsv,zone,evecsv,nstsv,h,nstsv,zzero,evecsvp, &
 209  nstsv)
 210 deallocate(igkignr,vgklnr,vgkcnr,gkcnr,tpgkcnr,vgqc,tpgqc,gqc,jlgqr)
 211 deallocate(evalsvp,evalsvnr,evecfv,evecsv,rwork)
 212 deallocate(h,vmat,apwalm,sfacgknr,ylmgq,sfacgq)
 213 deallocate(wfmt1,wfmt2,wfir1,wfir2)
 214 deallocate(zrhomt,zrhoir,zpchg,zvclmt,zvclir,zfmt,work,rfmt)
 215 return
 216 end subroutine
 217