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