1
2! Copyright (C) 2015 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 genvmmtsv(wfmt,vmat)
7use modmain
8use moddftu
9use modomp
10implicit none
11! arguments
12complex(8), intent(in) :: wfmt(npcmtmax,natmtot,nspinor,nstsv)
13complex(8), intent(inout) :: vmat(nstsv,nstsv)
14! local variables
15integer ist,jst,ispn,jspn
16integer is,ias,ld,nthd
17integer nrc,nrci,nrco,i
18integer nm,l,lm,npc,npci
19! allocatable arrays
20complex(8), allocatable :: wfmt1(:,:,:),wfmt2(:)
21! external functions
22complex(8), external :: zfmtinp
23if (.not.tvmatmt) return
24ld=lmmaxdm*nspinor
25allocate(wfmt1(npcmtmax,nspinor,nstsv),wfmt2(npcmtmax))
26! loop over atoms
27do ias=1,natmtot
28  is=idxis(ias)
29  nrc=nrcmt(is)
30  nrci=nrcmti(is)
31  nrco=nrc-nrci
32  npc=npcmt(is)
33  npci=npcmti(is)
34! convert wavefunctions to spherical harmonics
35  do ist=1,nstsv
36    do ispn=1,nspinor
37      call zfsht(nrc,nrci,wfmt(:,ias,ispn,ist),wfmt1(:,ispn,ist))
38    end do
39  end do
40! loop over second-variational states
41  do jst=1,nstsv
42! loop over spins
43    do ispn=1,nspinor
44      wfmt2(1:npc)=0.d0
45      do jspn=1,nspinor
46        do l=0,lmaxdm
47          if (tvmmt(l,ias)) then
48            nm=2*l+1
49            lm=l**2+1
50            if (l.le.lmaxi) then
51              call zgemm('N','N',nm,nrci,nm,zone,vmatmt(lm,ispn,lm,jspn,ias), &
52               ld,wfmt1(lm,jspn,jst),lmmaxi,zone,wfmt2(lm),lmmaxi)
53            end if
54            i=npci+lm
55            call zgemm('N','N',nm,nrco,nm,zone,vmatmt(lm,ispn,lm,jspn,ias),ld, &
56             wfmt1(i,jspn,jst),lmmaxo,zone,wfmt2(i),lmmaxo)
57          end if
58        end do
59      end do
60! compute the matrix elements
61      call holdthd(nstsv,nthd)
62!$OMP PARALLEL DO DEFAULT(SHARED) &
63!$OMP NUM_THREADS(nthd)
64      do ist=1,nstsv
65        vmat(ist,jst)=vmat(ist,jst)+zfmtinp(nrc,nrci,wrcmt(:,is), &
66         wfmt1(:,ispn,ist),wfmt2)
67      end do
68!$OMP END PARALLEL DO
69      call freethd(nthd)
70    end do
71  end do
72end do
73deallocate(wfmt1,wfmt2)
74end subroutine
75
76