1 2! Copyright (C) 2007-2008 J. K. Dewhurst, S. Sharma and E. K. U. Gross. 3! This file is distributed under the terms of the GNU General Public License. 4! See the file COPYING for license details. 5 6!BOP 7! !ROUTINE: genvclijjk 8! !INTERFACE: 9subroutine genvclijjk(ikp,vclijjk) 10! !USES: 11use modmain 12! !INPUT/OUTPUT PARAMETERS: 13! ikp : k-point from non-reduced set (in,integer) 14! vclijjk : Coulomb matrix elements (out,complex(nstsv,nstsv,nstsv,nkpt)) 15! !DESCRIPTION: 16! Calculates Coulomb matrix elements of the type $(i-jj-k)$. 17! 18! !REVISION HISTORY: 19! Created 2008 (Sharma) 20!EOP 21!BOC 22implicit none 23! arguments 24integer, intent(in) :: ikp 25complex(8), intent(out) :: vclijjk(nstsv,nstsv,nstsv,nkpt) 26! local variables 27integer ik,ist1,ist2,ist3 28integer iv(3),iq,ig 29real(8) vc(3) 30complex(8) z1 31! automatic arrays 32integer idx(nstsv) 33! allocatable arrays 34real(8), allocatable :: vgqc(:,:),gqc(:),gclgq(:),jlgqrmt(:,:,:) 35complex(8), allocatable :: apwalm(:,:,:,:),evecfv(:,:),evecsv(:,:) 36complex(8), allocatable :: ylmgq(:,:),sfacgq(:,:) 37complex(8), allocatable :: wfmt1(:,:,:,:),wfir1(:,:,:) 38complex(8), allocatable :: wfmt2(:,:,:,:),wfir2(:,:,:) 39complex(8), allocatable :: zrhomt(:,:,:),zrhoir(:,:) 40complex(8), allocatable :: zvclmt(:,:),zvclir(:) 41! external functions 42complex(8), external :: zfinp 43! allocate local arrays 44allocate(vgqc(3,ngvc),gqc(ngvc),gclgq(ngvc)) 45allocate(jlgqrmt(0:lnpsd,ngvc,nspecies)) 46allocate(apwalm(ngkmax,apwordmax,lmmaxapw,natmtot)) 47allocate(evecfv(nmatmax,nstfv),evecsv(nstsv,nstsv)) 48allocate(ylmgq(lmmaxo,ngvc),sfacgq(ngvc,natmtot)) 49allocate(wfmt1(npcmtmax,natmtot,nspinor,nstsv),wfir1(ngtc,nspinor,nstsv)) 50allocate(wfmt2(npcmtmax,natmtot,nspinor,nstsv),wfir2(ngtc,nspinor,nstsv)) 51allocate(zrhomt(npcmtmax,natmtot,nstsv),zrhoir(ngtc,nstsv)) 52allocate(zvclmt(npcmtmax,natmtot),zvclir(ngtc)) 53! get the eigenvectors from file for non-reduced k-point ikp 54call getevecfv(filext,0,vkl(:,ikp),vgkl(:,:,1,ikp),evecfv) 55call getevecsv(filext,0,vkl(:,ikp),evecsv) 56! find the matching coefficients 57call match(ngk(1,ikp),vgkc(:,:,1,ikp),gkc(:,1,ikp),sfacgk(:,:,1,ikp),apwalm) 58! index to all states 59do ist1=1,nstsv 60 idx(ist1)=ist1 61end do 62! calculate the wavefunctions for all states of passed non-reduced k-point ikp 63call genwfsv(.false.,.false.,nstsv,idx,ngdgc,igfc,ngk(1,ikp),igkig(:,1,ikp), & 64 apwalm,evecfv,evecsv,wfmt2,ngtc,wfir2) 65! start loop over reduced k-point set 66do ik=1,nkpt 67! determine the q-vector 68 iv(:)=ivk(:,ik)-ivk(:,ikp) 69 iv(:)=modulo(iv(:),ngridk(:)) 70! check if the q-point is in user-defined set 71 iv(:)=iv(:)*ngridq(:) 72 if (any(modulo(iv(:),ngridk(:)).ne.0)) cycle 73 iv(:)=iv(:)/ngridk(:) 74 iq=ivqiq(iv(1),iv(2),iv(3)) 75 vc(:)=vkc(:,ik)-vkc(:,ikp) 76 do ig=1,ngvc 77! determine the G+q-vectors 78 vgqc(:,ig)=vgc(:,ig)+vc(:) 79! G+q-vector length 80 gqc(ig)=sqrt(vgqc(1,ig)**2+vgqc(2,ig)**2+vgqc(3,ig)**2) 81! spherical harmonics for G+q-vectors 82 call genylmv(lmaxo,vgqc(:,ig),ylmgq(:,ig)) 83 end do 84! structure factors for G+q 85 call gensfacgp(ngvc,vgqc,ngvc,sfacgq) 86! generate the regularised Coulomb Green's function in G+q-space 87 call gengclgq(.true.,iq,ngvc,gqc,gclgq) 88! compute the required spherical Bessel functions 89 call genjlgprmt(lnpsd,ngvc,gqc,ngvc,jlgqrmt) 90! find the matching coefficients 91 call match(ngk(1,ik),vgkc(:,:,1,ik),gkc(:,1,ik),sfacgk(:,:,1,ik),apwalm) 92! get the eigenvectors from file 93 call getevecfv(filext,ik,vkl(:,ik),vgkl(:,:,:,ik),evecfv) 94 call getevecsv(filext,ik,vkl(:,ik),evecsv) 95! calculate the wavefunctions for all states of the reduced k-point 96 call genwfsv(.false.,.false.,nstsv,idx,ngdgc,igfc,ngk(1,ik),igkig(:,1,ik), & 97 apwalm,evecfv,evecsv,wfmt1,ngtc,wfir1) 98!----------------------------------------------! 99! valence-valence-valence contribution ! 100!----------------------------------------------! 101 do ist2=1,nstsv 102 do ist1=1,nstsv 103! calculate the complex overlap density for all states 104 call genzrho(.true.,.true.,ngtc,wfmt2(:,:,:,ist2),wfir2(:,:,ist2), & 105 wfmt1(:,:,:,ist1),wfir1(:,:,ist1),zrhomt(:,:,ist1),zrhoir(:,ist1)) 106 end do 107 do ist1=1,nstsv 108! compute the Coulomb potential 109 call genzvclmt(nrcmt,nrcmti,nrcmtmax,rlcmt,wprcmt,npcmtmax, & 110 zrhomt(:,:,ist1),zvclmt) 111 call zpotcoul(nrcmt,nrcmti,npcmt,npcmti,nrcmtmax,rlcmt,ngdgc,igfc,ngvc, & 112 gqc,gclgq,ngvc,jlgqrmt,ylmgq,sfacgq,zrhoir(:,ist1),npcmtmax,zvclmt, & 113 zvclir) 114 zvclir(:)=zvclir(:)*cfrc(:) 115 do ist3=ist1,nstsv 116 z1=zfinp(zrhomt(:,:,ist3),zrhoir(:,ist3),zvclmt,zvclir) 117 vclijjk(ist3,ist1,ist2,ik)=wqptnr*z1 118 end do 119 end do 120 end do 121! calculate the lower diagonal 122 do ist1=1,nstsv 123 do ist3=1,ist1-1 124 vclijjk(ist3,ist1,:,ik)=conjg(vclijjk(ist1,ist3,:,ik)) 125 end do 126 end do 127! end loop over reduced k-point set 128end do 129deallocate(vgqc,gqc,gclgq,jlgqrmt) 130deallocate(apwalm,evecfv,evecsv,ylmgq,sfacgq) 131deallocate(wfmt1,wfmt2,wfir1,wfir2) 132deallocate(zrhomt,zrhoir,zvclmt,zvclir) 133end subroutine 134!EOC 135 136