1
2! Copyright (C) 2014 K. Krieger, 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
6subroutine genvbmatk(vmt,vir,bmt,bir,ngp,igpig,wfmt,ld,wfgp,vbmat)
7use modmain
8use modomp
9implicit none
10! arguments
11! the potential and field are multiplied by the radial integration weights in
12! the muffin-tin and by the characteristic function in the interstitial region
13real(8), intent(in) :: vmt(npcmtmax,natmtot),vir(ngtot)
14real(8), intent(in) :: bmt(npcmtmax,natmtot,ndmag),bir(ngtot,ndmag)
15integer, intent(in) :: ngp(nspnfv),igpig(ngkmax,nspnfv)
16complex(8), intent(in) :: wfmt(npcmtmax,natmtot,nspinor,nstsv)
17integer, intent(in) :: ld
18complex(8), intent(in) :: wfgp(ld,nspinor,nstsv)
19complex(8), intent(out) :: vbmat(nstsv,nstsv)
20! local variables
21integer ist,jst,ispn,jspn
22integer is,ias,npc,igp,nthd
23! automatic arrays
24complex(8) wfmt2(npcmtmax,nspinor),z(ngkmax)
25! allocatable arrays
26complex(8), allocatable :: wfir1(:,:),wfir2(:,:)
27! external functions
28complex(8), external :: zdotc
29! zero the matrix elements
30vbmat(:,:)=0.d0
31!-------------------------!
32!     muffin-tin part     !
33!-------------------------!
34call holdthd(nstsv,nthd)
35!$OMP PARALLEL DO DEFAULT(SHARED) &
36!$OMP PRIVATE(wfmt2,ias,is,npc,ist) &
37!$OMP NUM_THREADS(nthd)
38do jst=1,nstsv
39  do ias=1,natmtot
40    is=idxis(ias)
41    npc=npcmt(is)
42! apply local potential and magnetic field to spinor wavefunction
43    if (ncmag) then
44! non-collinear case
45      call vbmk1(npc,vmt(:,ias),bmt(:,ias,1),bmt(:,ias,2),bmt(:,ias,3), &
46       wfmt(:,ias,1,jst),wfmt(:,ias,2,jst),wfmt2,wfmt2(:,2))
47    else
48! collinear case
49      call vbmk2(npc,vmt(:,ias),bmt(:,ias,1),wfmt(:,ias,1,jst), &
50       wfmt(:,ias,2,jst),wfmt2,wfmt2(:,2))
51    end if
52! compute the inner products
53    do ist=1,jst
54      vbmat(ist,jst)=vbmat(ist,jst) &
55       +zdotc(npc,wfmt(:,ias,1,ist),1,wfmt2,1) &
56       +zdotc(npc,wfmt(:,ias,2,ist),1,wfmt2(:,2),1)
57    end do
58  end do
59end do
60!$OMP END PARALLEL DO
61call freethd(nthd)
62!---------------------------!
63!     interstitial part     !
64!---------------------------!
65call holdthd(nstsv,nthd)
66!$OMP PARALLEL DEFAULT(SHARED) &
67!$OMP PRIVATE(wfir1,wfir2,z) &
68!$OMP PRIVATE(ispn,jspn,igp,ist) &
69!$OMP NUM_THREADS(nthd)
70allocate(wfir1(ngtot,nspinor),wfir2(ngtot,nspinor))
71!$OMP DO
72do jst=1,nstsv
73! Fourier transform wavefunction to real-space
74  do ispn=1,nspinor
75    jspn=jspnfv(ispn)
76    wfir1(:,ispn)=0.d0
77    do igp=1,ngp(jspn)
78      wfir1(igfft(igpig(igp,jspn)),ispn)=wfgp(igp,ispn,jst)
79    end do
80    call zfftifc(3,ngridg,1,wfir1(:,ispn))
81  end do
82! apply local potential and magnetic field to spinor wavefunction
83  if (ncmag) then
84! non-collinear case
85    call vbmk1(ngtot,vir,bir,bir(:,2),bir(:,3),wfir1,wfir1(:,2),wfir2, &
86     wfir2(:,2))
87  else
88! collinear case
89    call vbmk2(ngtot,vir,bir,wfir1,wfir1(:,2),wfir2,wfir2(:,2))
90  end if
91  do ispn=1,nspinor
92    jspn=jspnfv(ispn)
93! Fourier transform to G+p-space
94    call zfftifc(3,ngridg,-1,wfir2(:,ispn))
95    do igp=1,ngp(jspn)
96      z(igp)=wfir2(igfft(igpig(igp,jspn)),ispn)
97    end do
98    do ist=1,jst
99      vbmat(ist,jst)=vbmat(ist,jst)+zdotc(ngp(jspn),wfgp(:,ispn,ist),1,z,1)
100    end do
101  end do
102end do
103!$OMP END DO
104deallocate(wfir1,wfir2)
105!$OMP END PARALLEL
106call freethd(nthd)
107! lower triangular part
108do ist=1,nstsv
109  do jst=1,ist-1
110    vbmat(ist,jst)=conjg(vbmat(jst,ist))
111  end do
112end do
113return
114
115contains
116
117pure subroutine vbmk1(n,v,b1,b2,b3,wf11,wf12,wf21,wf22)
118implicit none
119! arguments
120integer, intent(in) :: n
121real(8), intent(in) :: v(n),b1(n),b2(n),b3(n)
122complex(8), intent(in) :: wf11(n),wf12(n)
123complex(8), intent(out) :: wf21(n),wf22(n)
124! local variables
125integer i
126complex(8) z1
127do i=1,n
128  z1=cmplx(b1(i),b2(i),8)
129  wf21(i)=(v(i)+b3(i))*wf11(i)+conjg(z1)*wf12(i)
130  wf22(i)=(v(i)-b3(i))*wf12(i)+z1*wf11(i)
131end do
132end subroutine
133
134pure subroutine vbmk2(n,v,b,wf11,wf12,wf21,wf22)
135implicit none
136! arguments
137integer, intent(in) :: n
138real(8), intent(in) :: v(n),b(n)
139complex(8), intent(in) :: wf11(n),wf12(n)
140complex(8), intent(out) :: wf21(n),wf22(n)
141! local variables
142integer i
143do i=1,n
144  wf21(i)=(v(i)+b(i))*wf11(i)
145  wf22(i)=(v(i)-b(i))*wf12(i)
146end do
147end subroutine
148
149end subroutine
150
151