File: dhmlrad.f90

package info (click to toggle)
elkcode 5.4.24-2
  • links: PTS, VCS
  • area: main
  • in suites: buster, sid
  • size: 12,840 kB
  • sloc: f90: 48,415; fortran: 22,457; perl: 965; makefile: 384; sh: 369; python: 105; ansic: 67
file content (174 lines) | stat: -rw-r--r-- 5,084 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

! Copyright (C) 2013 J. K. Dewhurst, S. Sharma and E. K. U. Gross.
! This file is distributed under the terms of the GNU General Public License.
! See the file COPYING for license details.

subroutine dhmlrad
use modmain
use modphonon
implicit none
! local variables
integer is,ias
integer nr,nri,nro,iro,ir
integer l1,l2,l3,m2,lm2
integer npi,i,ilo,jlo,io,jo
real(8) t1,t2
! automatic arrays
real(8) fr1(nrmtmax),fr2(nrmtmax)
! external functions
real(8) splint
external splint
! begin loops over atoms and species
do ias=1,natmtot
  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)
          lm2=0
          do l2=0,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)
                fr1(ir)=t1*dble(dvsmt(i,ias))
                fr2(ir)=t1*aimag(dvsmt(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)
                fr1(ir)=t1*dble(dvsmt(i,ias))
                fr2(ir)=t1*aimag(dvsmt(i,ias))
                i=i+lmmaxo
              end do
              t1=splint(nr,rsp(:,is),fr1)
              t2=splint(nr,rsp(:,is),fr2)
              dhaa(lm2,jo,l3,io,l1,ias)=cmplx(t1,t2,8)
            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)
                fr1(ir)=t1*dble(dvsmt(i,ias))
                fr2(ir)=t1*aimag(dvsmt(i,ias))
                i=i+lmmaxo
              end do
              t1=splint(nro,rsp(iro,is),fr1(iro))
              t2=splint(nro,rsp(iro,is),fr2(iro))
              dhaa(lm2,jo,l3,io,l1,ias)=cmplx(t1,t2,8)
            end do
          end do
        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)
        lm2=0
        do l2=0,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)
              fr1(ir)=t1*dble(dvsmt(i,ias))
              fr2(ir)=t1*aimag(dvsmt(i,ias))
              i=i+lmmaxi
            end do
            do ir=iro,nr
              t1=lofr(ir,1,ilo,ias)*apwfr(ir,1,io,l3,ias)*r2sp(ir,is)
              fr1(ir)=t1*dble(dvsmt(i,ias))
              fr2(ir)=t1*aimag(dvsmt(i,ias))
              i=i+lmmaxo
            end do
            t1=splint(nr,rsp(:,is),fr1)
            t2=splint(nr,rsp(:,is),fr2)
            dhloa(lm2,io,l3,ilo,ias)=cmplx(t1,t2,8)
          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)
              fr1(ir)=t1*dble(dvsmt(i,ias))
              fr2(ir)=t1*aimag(dvsmt(i,ias))
              i=i+lmmaxo
            end do
            t1=splint(nro,rsp(iro,is),fr1(iro))
            t2=splint(nro,rsp(iro,is),fr2(iro))
            dhloa(lm2,io,l3,ilo,ias)=cmplx(t1,t2,8)
          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)
      lm2=0
      do l2=0,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)
            fr1(ir)=t1*dble(dvsmt(i,ias))
            fr2(ir)=t1*aimag(dvsmt(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)
            fr1(ir)=t1*dble(dvsmt(i,ias))
            fr2(ir)=t1*aimag(dvsmt(i,ias))
            i=i+lmmaxo
          end do
          t1=splint(nr,rsp(:,is),fr1)
          t2=splint(nr,rsp(:,is),fr2)
          dhlolo(lm2,jlo,ilo,ias)=cmplx(t1,t2,8)
        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)
            fr1(ir)=t1*dble(dvsmt(i,ias))
            fr2(ir)=t1*aimag(dvsmt(i,ias))
            i=i+lmmaxo
          end do
          t1=splint(nro,rsp(iro,is),fr1(iro))
          t2=splint(nro,rsp(iro,is),fr2(iro))
          dhlolo(lm2,jlo,ilo,ias)=cmplx(t1,t2,8)
        end do
      end do
    end do
  end do
! end loops over atoms and species
end do
return
end subroutine