1
2! Copyright (C) 2002-2016 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
6!BOP
7! !ROUTINE: hmlrad
8! !INTERFACE:
9subroutine hmlrad
10! !USES:
11use modmain
12use modomp
13! !DESCRIPTION:
14!   Calculates the radial Hamiltonian integrals of the APW and local-orbital
15!   basis functions. In other words, for atom $\alpha$, it computes integrals of
16!   the form
17!   $$ h^{\alpha}_{qq';ll'l''m''}=\begin{cases}
18!    \int_0^{R_i}u^{\alpha}_{q;l}(r)H u^{\alpha}_{q';l'}(r)r^2dr & l''=0 \\
19!    \int_0^{R_i}u^{\alpha}_{q;l}(r)V^{\alpha}_{l''m''}(r)
20!    u^{\alpha}_{q';l'}(r)r^2dr & l''>0 \end{cases}, $$
21!   where $u^{\alpha}_{q;l}$ is the $q$th APW radial function for angular
22!   momentum $l$; $H$ is the Hamiltonian of the radial Schr\"{o}dinger equation;
23!   and $V^{\alpha}_{l''m''}$ is the muffin-tin Kohn-Sham potential. Similar
24!   integrals are calculated for APW-local-orbital and
25!   local-orbital-local-orbital contributions.
26!
27! !REVISION HISTORY:
28!   Created December 2003 (JKD)
29!   Updated for compressed muffin-tin functions, March 2016 (JKD)
30!EOP
31!BOC
32implicit none
33! local variables
34integer is,ias,nthd
35integer nr,nri,iro
36integer ir,npi,i
37integer l1,l2,l3,m2,lm2
38integer io,jo,ilo,jlo
39real(8) sm,t1
40! begin loops over atoms and species
41call holdthd(natmtot,nthd)
42!$OMP PARALLEL DO DEFAULT(SHARED) &
43!$OMP PRIVATE(is,nr,nri,iro,npi) &
44!$OMP PRIVATE(l1,l2,l3,io,jo,sm) &
45!$OMP PRIVATE(lm2,m2,i,ir,t1,ilo,jlo) &
46!$OMP NUM_THREADS(nthd)
47do ias=1,natmtot
48  is=idxis(ias)
49  nr=nrmt(is)
50  nri=nrmti(is)
51  iro=nri+1
52  npi=npmti(is)
53!---------------------------!
54!     APW-APW integrals     !
55!---------------------------!
56  do l1=0,lmaxapw
57    do io=1,apword(l1,is)
58      do l3=0,lmaxapw
59        do jo=1,apword(l3,is)
60          if (l1.eq.l3) then
61            sm=sum(apwfr(1:nr,1,io,l1,ias)*apwfr(1:nr,2,jo,l3,ias) &
62             *wrmt(1:nr,is))
63            haa(1,jo,l3,io,l1,ias)=sm/y00
64          else
65            haa(1,jo,l3,io,l1,ias)=0.d0
66          end if
67          if (l1.ge.l3) then
68            lm2=1
69            do l2=1,lmaxi
70              do m2=-l2,l2
71                lm2=lm2+1
72                sm=0.d0
73                i=lm2
74                do ir=1,nri
75                  t1=apwfr(ir,1,io,l1,ias)*apwfr(ir,1,jo,l3,ias)*wrmt(ir,is)
76                  sm=sm+t1*vsmt(i,ias)
77                  i=i+lmmaxi
78                end do
79                do ir=iro,nr
80                  t1=apwfr(ir,1,io,l1,ias)*apwfr(ir,1,jo,l3,ias)*wrmt(ir,is)
81                  sm=sm+t1*vsmt(i,ias)
82                  i=i+lmmaxo
83                end do
84                haa(lm2,jo,l3,io,l1,ias)=sm
85                haa(lm2,io,l1,jo,l3,ias)=sm
86              end do
87            end do
88            do l2=lmaxi+1,lmaxo
89              do m2=-l2,l2
90                lm2=lm2+1
91                sm=0.d0
92                i=npi+lm2
93                do ir=iro,nr
94                  t1=apwfr(ir,1,io,l1,ias)*apwfr(ir,1,jo,l3,ias)*wrmt(ir,is)
95                  sm=sm+t1*vsmt(i,ias)
96                  i=i+lmmaxo
97                end do
98                haa(lm2,jo,l3,io,l1,ias)=sm
99                haa(lm2,io,l1,jo,l3,ias)=sm
100              end do
101            end do
102          end if
103        end do
104      end do
105    end do
106  end do
107!-------------------------------------!
108!     local-orbital-APW integrals     !
109!-------------------------------------!
110  do ilo=1,nlorb(is)
111    l1=lorbl(ilo,is)
112    do l3=0,lmaxapw
113      do io=1,apword(l3,is)
114        if (l1.eq.l3) then
115          sm=sum(lofr(1:nr,1,ilo,ias)*apwfr(1:nr,2,io,l3,ias)*wrmt(1:nr,is))
116          hloa(1,io,l3,ilo,ias)=sm/y00
117        else
118          hloa(1,io,l3,ilo,ias)=0.d0
119        end if
120        lm2=1
121        do l2=1,lmaxi
122          do m2=-l2,l2
123            lm2=lm2+1
124            sm=0.d0
125            i=lm2
126            do ir=1,nri
127              t1=lofr(ir,1,ilo,ias)*apwfr(ir,1,io,l3,ias)*wrmt(ir,is)
128              sm=sm+t1*vsmt(i,ias)
129              i=i+lmmaxi
130            end do
131            do ir=nri+1,nr
132              t1=lofr(ir,1,ilo,ias)*apwfr(ir,1,io,l3,ias)*wrmt(ir,is)
133              sm=sm+t1*vsmt(i,ias)
134              i=i+lmmaxo
135            end do
136            hloa(lm2,io,l3,ilo,ias)=sm
137          end do
138        end do
139        do l2=lmaxi+1,lmaxo
140          do m2=-l2,l2
141            lm2=lm2+1
142            sm=0.d0
143            i=npi+lm2
144            do ir=iro,nr
145              t1=lofr(ir,1,ilo,ias)*apwfr(ir,1,io,l3,ias)*wrmt(ir,is)
146              sm=sm+t1*vsmt(i,ias)
147              i=i+lmmaxo
148            end do
149            hloa(lm2,io,l3,ilo,ias)=sm
150          end do
151        end do
152      end do
153    end do
154  end do
155!-----------------------------------------------!
156!     local-orbital-local-orbital integrals     !
157!-----------------------------------------------!
158  do ilo=1,nlorb(is)
159    l1=lorbl(ilo,is)
160    do jlo=1,nlorb(is)
161      l3=lorbl(jlo,is)
162      if (l1.eq.l3) then
163        sm=sum(lofr(1:nr,1,ilo,ias)*lofr(1:nr,2,jlo,ias)*wrmt(1:nr,is))
164        hlolo(1,jlo,ilo,ias)=sm/y00
165      else
166        hlolo(1,jlo,ilo,ias)=0.d0
167      end if
168      lm2=1
169      do l2=1,lmaxi
170        do m2=-l2,l2
171          lm2=lm2+1
172          sm=0.d0
173          i=lm2
174          do ir=1,nri
175            t1=lofr(ir,1,ilo,ias)*lofr(ir,1,jlo,ias)*wrmt(ir,is)
176            sm=sm+t1*vsmt(i,ias)
177            i=i+lmmaxi
178          end do
179          do ir=iro,nr
180            t1=lofr(ir,1,ilo,ias)*lofr(ir,1,jlo,ias)*wrmt(ir,is)
181            sm=sm+t1*vsmt(i,ias)
182            i=i+lmmaxo
183          end do
184          hlolo(lm2,jlo,ilo,ias)=sm
185        end do
186      end do
187      do l2=lmaxi+1,lmaxo
188        do m2=-l2,l2
189          lm2=lm2+1
190          sm=0.d0
191          i=npi+lm2
192          do ir=iro,nr
193            t1=lofr(ir,1,ilo,ias)*lofr(ir,1,jlo,ias)*wrmt(ir,is)
194            sm=sm+t1*vsmt(i,ias)
195            i=i+lmmaxo
196          end do
197          hlolo(lm2,jlo,ilo,ias)=sm
198        end do
199      end do
200    end do
201  end do
202! end loops over atoms and species
203end do
204!$OMP END PARALLEL DO
205call freethd(nthd)
206end subroutine
207!EOC
208
209