1
2! Copyright (C) 2007 J. K. Dewhurst, S. Sharma and C. Ambrosch-Draxl.
3! This file is distributed under the terms of the GNU General Public License.
4! See the file COPYING for license details.
5
6subroutine genvmatk(vmt,vir,ngp,igpig,wfmt,ld,wfgp,vmat)
7use modmain
8use modomp
9implicit none
10! arguments
11! the potential is multiplied by the radial integration weights in the
12! muffin-tin and by the characteristic function in the interstitial region
13real(8), intent(in) :: vmt(npcmtmax,natmtot),vir(ngtot)
14integer, intent(in) :: ngp(nspnfv),igpig(ngkmax,nspnfv)
15complex(8), intent(in) :: wfmt(npcmtmax,natmtot,nspinor,nstsv)
16integer, intent(in) :: ld
17complex(8), intent(in) :: wfgp(ld,nspinor,nstsv)
18complex(8), intent(out) :: vmat(nstsv,nstsv)
19! local variables
20integer ist,jst,ispn,jspn
21integer is,ias,npc,igp,nthd
22! automatic arrays
23complex(8) wfmt1(npcmtmax),z(ngkmax)
24! allocatable arrays
25complex(8), allocatable :: wfir(:)
26! external functions
27complex(8), external :: zdotc
28! zero the matrix elements
29vmat(:,:)=0.d0
30!-------------------------!
31!     muffin-tin part     !
32!-------------------------!
33call holdthd(nstsv,nthd)
34!$OMP PARALLEL DO DEFAULT(SHARED) &
35!$OMP PRIVATE(wfmt1,ispn,ias) &
36!$OMP PRIVATE(is,npc,ist) &
37!$OMP NUM_THREADS(nthd)
38do jst=1,nstsv
39  do ispn=1,nspinor
40    do ias=1,natmtot
41      is=idxis(ias)
42      npc=npcmt(is)
43! apply potential to wavefunction
44      wfmt1(1:npc)=vmt(1:npc,ias)*wfmt(1:npc,ias,ispn,jst)
45! compute the inner products
46      do ist=1,jst
47        vmat(ist,jst)=vmat(ist,jst)+zdotc(npc,wfmt(:,ias,ispn,ist),1,wfmt1,1)
48      end do
49    end do
50  end do
51end do
52!$OMP END PARALLEL DO
53call freethd(nthd)
54!---------------------------!
55!     interstitial part     !
56!---------------------------!
57call holdthd(nstsv,nthd)
58!$OMP PARALLEL DEFAULT(SHARED) &
59!$OMP PRIVATE(wfir,z,ispn,jspn) &
60!$OMP PRIVATE(igp,ist) &
61!$OMP NUM_THREADS(nthd)
62allocate(wfir(ngtot))
63!$OMP DO
64do jst=1,nstsv
65  do ispn=1,nspinor
66    jspn=jspnfv(ispn)
67! Fourier transform wavefunction to real-space
68    wfir(:)=0.d0
69    do igp=1,ngp(jspn)
70      wfir(igfft(igpig(igp,jspn)))=wfgp(igp,ispn,jst)
71    end do
72    call zfftifc(3,ngridg,1,wfir)
73! apply potential to wavefunction
74    wfir(1:ngtot)=vir(1:ngtot)*wfir(1:ngtot)
75! Fourier transform to G+p-space
76    call zfftifc(3,ngridg,-1,wfir)
77    do igp=1,ngp(jspn)
78      z(igp)=wfir(igfft(igpig(igp,jspn)))
79    end do
80    do ist=1,jst
81! compute inner product
82      vmat(ist,jst)=vmat(ist,jst)+zdotc(ngp(jspn),wfgp(:,ispn,ist),1,z,1)
83    end do
84  end do
85end do
86!$OMP END DO
87deallocate(wfir)
88!$OMP END PARALLEL
89call freethd(nthd)
90! lower triangular part
91do ist=1,nstsv
92  do jst=1,ist-1
93    vmat(ist,jst)=conjg(vmat(jst,ist))
94  end do
95end do
96end subroutine
97
98