File: hmlrad.f90

package info (click to toggle)
elkcode 5.4.24-2
  • links: PTS, VCS
  • area: main
  • in suites: buster
  • size: 12,840 kB
  • sloc: f90: 48,415; fortran: 22,457; perl: 965; makefile: 384; sh: 369; python: 105; ansic: 67
file content (227 lines) | stat: -rw-r--r-- 6,549 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227

! Copyright (C) 2002-2016 J. K. Dewhurst, S. Sharma and C. Ambrosch-Draxl.
! This file is distributed under the terms of the GNU General Public License.
! See the file COPYING for license details.

!BOP
! !ROUTINE: hmlrad
! !INTERFACE:
subroutine hmlrad
! !USES:
use modmain
use modomp
! !DESCRIPTION:
!   Calculates the radial Hamiltonian integrals of the APW and local-orbital
!   basis functions. In other words, for atom $\alpha$, it computes integrals of
!   the form
!   $$ h^{\alpha}_{qq';ll'l''m''}=\begin{cases}
!    \int_0^{R_i}u^{\alpha}_{q;l}(r)H u^{\alpha}_{q';l'}(r)r^2dr & l''=0 \\
!    \int_0^{R_i}u^{\alpha}_{q;l}(r)V^{\alpha}_{l''m''}(r)
!    u^{\alpha}_{q';l'}(r)r^2dr & l''>0 \end{cases}, $$
!   where $u^{\alpha}_{q;l}$ is the $q$th APW radial function for angular
!   momentum $l$; $H$ is the Hamiltonian of the radial Schr\"{o}dinger equation;
!   and $V^{\alpha}_{l''m''}$ is the muffin-tin Kohn-Sham potential. Similar
!   integrals are calculated for APW-local-orbital and
!   local-orbital-local-orbital contributions.
!
! !REVISION HISTORY:
!   Created December 2003 (JKD)
!   Updated for compressed muffin-tin functions, March 2016 (JKD)
!EOP
!BOC
implicit none
! local variables
integer is,ias,nthd
integer nr,nri,nro
integer iro,ir,npi,i
integer l1,l2,l3,m2,lm2
integer io,jo,ilo,jlo
real(8) t1
! allocatable arrays
real(8), allocatable :: fr(:)
! external functions
real(8) splint
external splint
! begin loops over atoms and species
call omp_hold(natmtot,nthd)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(fr,is,nr,nri,nro,iro,npi) &
!$OMP PRIVATE(l1,l2,l3,io,jo,ir,t1) &
!$OMP PRIVATE(lm2,m2,i,ilo,jlo) &
!$OMP NUM_THREADS(nthd)
!$OMP DO
do ias=1,natmtot
  allocate(fr(nrmtmax))
  is=idxis(ias)
  nr=nrmt(is)
  nri=nrmti(is)
  nro=nr-nri
  iro=nri+1
  npi=npmti(is)
!---------------------------!
!     APW-APW integrals     !
!---------------------------!
  do l1=0,lmaxapw
    do io=1,apword(l1,is)
      do l3=0,lmaxapw
        do jo=1,apword(l3,is)
          if (l1.eq.l3) then
            do ir=1,nr
              fr(ir)=apwfr(ir,1,io,l1,ias)*apwfr(ir,2,jo,l3,ias)*r2sp(ir,is)
            end do
            t1=splint(nr,rsp(:,is),fr)
            haa(1,jo,l3,io,l1,ias)=t1/y00
          else
            haa(1,jo,l3,io,l1,ias)=0.d0
          end if
          if (l1.ge.l3) then
            lm2=1
            do l2=1,lmaxi
              do m2=-l2,l2
                lm2=lm2+1
                i=lm2
                do ir=1,nri
                  t1=apwfr(ir,1,io,l1,ias)*apwfr(ir,1,jo,l3,ias)*r2sp(ir,is)
                  fr(ir)=t1*vsmt(i,ias)
                  i=i+lmmaxi
                end do
                do ir=iro,nr
                  t1=apwfr(ir,1,io,l1,ias)*apwfr(ir,1,jo,l3,ias)*r2sp(ir,is)
                  fr(ir)=t1*vsmt(i,ias)
                  i=i+lmmaxo
                end do
                t1=splint(nr,rsp(:,is),fr)
                haa(lm2,jo,l3,io,l1,ias)=t1
                haa(lm2,io,l1,jo,l3,ias)=t1
              end do
            end do
            do l2=lmaxi+1,lmaxo
              do m2=-l2,l2
                lm2=lm2+1
                i=npi+lm2
                do ir=iro,nr
                  t1=apwfr(ir,1,io,l1,ias)*apwfr(ir,1,jo,l3,ias)*r2sp(ir,is)
                  fr(ir)=t1*vsmt(i,ias)
                  i=i+lmmaxo
                end do
                t1=splint(nro,rsp(iro,is),fr(iro))
                haa(lm2,jo,l3,io,l1,ias)=t1
                haa(lm2,io,l1,jo,l3,ias)=t1
              end do
            end do
          end if
        end do
      end do
    end do
  end do
!-------------------------------------!
!     local-orbital-APW integrals     !
!-------------------------------------!
  do ilo=1,nlorb(is)
    l1=lorbl(ilo,is)
    do l3=0,lmaxapw
      do io=1,apword(l3,is)
        if (l1.eq.l3) then
          do ir=1,nr
            fr(ir)=lofr(ir,1,ilo,ias)*apwfr(ir,2,io,l3,ias)*r2sp(ir,is)
          end do
          t1=splint(nr,rsp(:,is),fr)
          hloa(1,io,l3,ilo,ias)=t1/y00
        else
          hloa(1,io,l3,ilo,ias)=0.d0
        end if
        lm2=1
        do l2=1,lmaxi
          do m2=-l2,l2
            lm2=lm2+1
            i=lm2
            do ir=1,nri
              t1=lofr(ir,1,ilo,ias)*apwfr(ir,1,io,l3,ias)*r2sp(ir,is)
              fr(ir)=t1*vsmt(i,ias)
              i=i+lmmaxi
            end do
            do ir=nri+1,nr
              t1=lofr(ir,1,ilo,ias)*apwfr(ir,1,io,l3,ias)*r2sp(ir,is)
              fr(ir)=t1*vsmt(i,ias)
              i=i+lmmaxo
            end do
            t1=splint(nr,rsp(:,is),fr)
            hloa(lm2,io,l3,ilo,ias)=t1
          end do
        end do
        do l2=lmaxi+1,lmaxo
          do m2=-l2,l2
            lm2=lm2+1
            i=npi+lm2
            do ir=iro,nr
              t1=lofr(ir,1,ilo,ias)*apwfr(ir,1,io,l3,ias)*r2sp(ir,is)
              fr(ir)=t1*vsmt(i,ias)
              i=i+lmmaxo
            end do
            t1=splint(nro,rsp(iro,is),fr(iro))
            hloa(lm2,io,l3,ilo,ias)=t1
          end do
        end do
      end do
    end do
  end do
!-----------------------------------------------!
!     local-orbital-local-orbital integrals     !
!-----------------------------------------------!
  do ilo=1,nlorb(is)
    l1=lorbl(ilo,is)
    do jlo=1,nlorb(is)
      l3=lorbl(jlo,is)
      if (l1.eq.l3) then
        do ir=1,nr
          fr(ir)=lofr(ir,1,ilo,ias)*lofr(ir,2,jlo,ias)*r2sp(ir,is)
        end do
        t1=splint(nr,rsp(:,is),fr)
        hlolo(1,jlo,ilo,ias)=t1/y00
      else
        hlolo(1,jlo,ilo,ias)=0.d0
      end if
      lm2=1
      do l2=1,lmaxi
        do m2=-l2,l2
          lm2=lm2+1
          i=lm2
          do ir=1,nri
            t1=lofr(ir,1,ilo,ias)*lofr(ir,1,jlo,ias)*r2sp(ir,is)
            fr(ir)=t1*vsmt(i,ias)
            i=i+lmmaxi
          end do
          do ir=iro,nr
            t1=lofr(ir,1,ilo,ias)*lofr(ir,1,jlo,ias)*r2sp(ir,is)
            fr(ir)=t1*vsmt(i,ias)
            i=i+lmmaxo
          end do
          t1=splint(nr,rsp(:,is),fr)
          hlolo(lm2,jlo,ilo,ias)=t1
        end do
      end do
      do l2=lmaxi+1,lmaxo
        do m2=-l2,l2
          lm2=lm2+1
          i=npi+lm2
          do ir=iro,nr
            t1=lofr(ir,1,ilo,ias)*lofr(ir,1,jlo,ias)*r2sp(ir,is)
            fr(ir)=t1*vsmt(i,ias)
            i=i+lmmaxo
          end do
          t1=splint(nro,rsp(iro,is),fr(iro))
          hlolo(lm2,jlo,ilo,ias)=t1
        end do
      end do
    end do
  end do
  deallocate(fr)
! end loops over atoms and species
end do
!$OMP END DO
!$OMP END PARALLEL
call omp_free(nthd)
return
end subroutine
!EOC