Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 : MODULE m_mkgylm
7 : CONTAINS
8 206921 : SUBROUTINE mkgylm(jspins,rv,thet,nsp, rh,rhdr,rhdt,rhdf,rhdrr,rhdtt,rhdff, rhdtf,rhdrt,rhdrf,grad,kt)
9 : !c.....------------------------------------------------------------------
10 : !c by use of charge density and its polar coord. gradient components
11 : !c c calculate agr and others used to evaluate gradient
12 : !c c contributions to potential and energy. t.a. 1996.
13 : !c.....------------------------------------------------------------------
14 : !c ro=sum(ro*ylh), rdr=sum(drr*ylh), drdr=sum(ddrr*ylh),
15 : !c c rdt=sum(ro*ylht1), rdtt=sum(ro*ylht2), ...
16 : !c c rdf=sum(ro*ylhf1), rdff=sum(ro*ylhf2), ...
17 : !c c rdtf=sum(ro*ylhtf), rdff=sum(ro*ylhf2), ...
18 : !c c rdrt=sum(drr*ylht1),rdrf=sum(drr*ylhf1),!
19 : !c agr: abs(grad(ro)), g2r: laplacian(ro),
20 : !c c gggr: grad(ro)*grad(agr),
21 : !c c grgru,d: grad(ro)*grad(rou),for rod., gzgr: grad(zeta)*grad(ro).
22 : !
23 : !c dagrr,-t,-f: d(agr)/dr, d(agr)/dth/r, d(agr)/dfi/r/sint.
24 : !c.....------------------------------------------------------------------
25 : USE m_types
26 : IMPLICIT NONE
27 : !C ..
28 : !C .. Arguments ..
29 : REAL, INTENT (IN) :: rv
30 : REAL, INTENT (IN) :: thet(:)
31 : REAL, INTENT (IN) :: rh(:,:),rhdr(:,:)
32 : REAL, INTENT (IN) :: rhdf(:,:),rhdrr(:,:)
33 : REAL, INTENT (IN) :: rhdtt(:,:),rhdff(:,:)
34 : REAL, INTENT (IN) :: rhdtf(:,:),rhdrt(:,:)
35 : REAL, INTENT (IN) :: rhdrf(:,:),rhdt(:,:)
36 : TYPE(t_gradients),INTENT(INOUT) :: grad
37 : INTEGER,INTENT(IN) :: jspins,nsp
38 : INTEGER,INTENT(IN) :: kt !index of first point to use in gradients
39 :
40 : !C ..
41 : !C .. Locals ..
42 : INTEGER i,js,fac
43 : REAL chsml, dagrf,dagrfd,dagrfu,dagrr,dagrrd,dagrru,dagrt,&
44 : dagrtd,dagrtu,drdr,dzdfs,dzdr,dzdtr,grf,grfd,grfu,grr,&
45 : grrd,grru,grt,grtd,grtu,rdf,rdfd,rdff,rdffd,rdffu,rdfu,&
46 : rdr,rdrd,rdrf,rdrfd,rdrfu,rdrrd,rdrru,rdrt,rdrtd,rdrtu,&
47 : rdru,rdt,rdtd,rdtf,rdtfd,rdtfu,rdtt,rdttd,rdttu,rdtu,&
48 : ro,ro2,rod,rou,rv1,rv2,rv3,rvsin1,sint1,sint2,tant1
49 :
50 206921 : chsml = 1.e-10
51 :
52 206921 : rv1 = rv
53 206921 : rv2 = rv1**2
54 206921 : rv3 = rv1**3
55 :
56 206921 : IF (ALLOCATED(grad%gr)) THEN
57 : ! Gradients for libxc
58 158200 : DO js=1,jspins
59 19368200 : DO i=1,nsp
60 19323000 : grad%gr(:,kt+i,js)=[rhdr(i,js),rhdt(i,js),rhdf(i,js)]
61 : ENDDO
62 : ENDDO
63 : !contracted gradients for libxc
64 45200 : IF (ALLOCATED(grad%sigma)) THEN
65 : ! Contracted gradients for libxc
66 22600 : IF (jspins==1) THEN
67 0 : DO i=1,nsp
68 0 : grad%sigma(1,kt+i)= dot_PRODUCT(grad%gr(:,kt+i,1),grad%gr(:,kt+i,1))
69 : ENDDO
70 : ELSE
71 3864600 : DO i=1,nsp
72 3842000 : grad%sigma(1,kt+i)= dot_PRODUCT(grad%gr(:,kt+i,1),grad%gr(:,kt+i,1))
73 3842000 : grad%sigma(2,kt+i)= dot_PRODUCT(grad%gr(:,kt+i,1),grad%gr(:,kt+i,2))
74 15390600 : grad%sigma(3,kt+i)= dot_PRODUCT(grad%gr(:,kt+i,2),grad%gr(:,kt+i,2))
75 : ENDDO
76 : ENDIF
77 : END IF
78 45200 : IF (ALLOCATED(grad%laplace)) THEN
79 : !Lapacian of density
80 : !fac=MERGE(2,1,jspins==1)
81 : fac=1
82 113000 : DO js=1,jspins
83 7751800 : DO i=1,nsp
84 : grad%laplace(kt+i,js) = (rhdrr(i,js) + 2.e0*rhdr(i,js)/rv1 +&
85 7729200 : (rhdtt(i,js)+rhdt(i,js)/TAN(thet(i))+rhdff(i,js)/SIN(thet(i))**2)/rv2)/fac
86 : ENDDO
87 : ENDDO
88 : ENDIF
89 : RETURN !Do not calculate arrays for in-build GGA
90 : END IF
91 : ! Old code for in-build xcpots
92 161721 : IF(allocated(grad%agrt)) THEN
93 161721 : IF (jspins==1) THEN
94 :
95 7766442 : points_1 : DO i = 1,nsp
96 :
97 7729260 : grad%agrt(kt+i) = 0.0
98 7729260 : grad%agru(kt+i) = 0.0
99 7729260 : grad%agrd(kt+i) = 0.0
100 7729260 : grad%g2rt(kt+i) = 0.0
101 7729260 : grad%g2ru(kt+i) = 0.0
102 7729260 : grad%g2rd(kt+i) = 0.0
103 7729260 : grad%gggrt(kt+i) = 0.0
104 7729260 : grad%gggru(kt+i) = 0.0
105 7729260 : grad%gggrd(kt+i) = 0.0
106 7729260 : grad%grgru(kt+i) = 0.0
107 7729260 : grad%grgrd(kt+i) = 0.0
108 7729260 : grad%gzgr(kt+i) = 0.0
109 :
110 7729260 : ro = rh(i,1)
111 :
112 7729260 : IF (ro<chsml) CYCLE points_1
113 7729260 : sint1 = sin(thet(i))
114 7729260 : sint2 = sint1**2
115 7729260 : tant1 = tan(thet(i))
116 7729260 : rvsin1 = rv1*sint1
117 :
118 7729260 : rou = ro/2
119 7729260 : rdru = rhdr(i,1)/2
120 7729260 : rdtu = rhdt(i,1)/2
121 7729260 : rdfu = rhdf(i,1)/2
122 7729260 : rdrru = rhdrr(i,1)/2
123 7729260 : rdttu = rhdtt(i,1)/2
124 7729260 : rdffu = rhdff(i,1)/2
125 7729260 : rdtfu = rhdtf(i,1)/2
126 7729260 : rdrtu = rhdrt(i,1)/2
127 7729260 : rdrfu = rhdrf(i,1)/2
128 :
129 7729260 : rod = rou
130 7729260 : rdrd = rdru
131 7729260 : rdtd = rdtu
132 7729260 : rdfd = rdfu
133 7729260 : rdrrd = rdrru
134 7729260 : rdttd = rdttu
135 7729260 : rdffd = rdffu
136 7729260 : rdtfd = rdtfu
137 7729260 : rdrtd = rdrtu
138 7729260 : rdrfd = rdrfu
139 :
140 7729260 : rdr = rdru + rdrd
141 7729260 : rdt = rdtu + rdtd
142 7729260 : rdf = rdfu + rdfd
143 7729260 : drdr = rdrru + rdrrd
144 7729260 : rdtt = rdttu + rdttd
145 7729260 : rdff = rdffu + rdffd
146 7729260 : rdtf = rdtfu + rdtfd
147 7729260 : rdrt = rdrtu + rdrtd
148 7729260 : rdrf = rdrfu + rdrfd
149 :
150 7729260 : ro2 = ro**2
151 :
152 7729260 : grr = rdr
153 7729260 : grt = rdt/rv1
154 7729260 : grf = rdf/rvsin1
155 :
156 7729260 : grad%agrt(kt+i) = sqrt(grr**2+grt**2+grf**2)
157 :
158 7729260 : IF (grad%agrt(kt+i)<chsml) CYCLE points_1
159 :
160 : dagrr = (rdr*drdr*rv3+rdt* (rdrt*rv1-rdt)+ &
161 7729260 : rdf* (rdrf*rv1-rdf)/sint2)/grad%agrt(kt+i)/rv3
162 :
163 : dagrt =(rdr*rdrt*rv2+rdt*rdtt+rdf* (-rdf/tant1+rdtf)/sint2)/&
164 7729260 : (grad%agrt(kt+i)*rv3)
165 :
166 : dagrf = (rdr*rdrf*rv2+rdt*rdtf+rdf*rdff/sint2)/&
167 7729260 : (grad%agrt(kt+i)*rv3*sint1)
168 :
169 : grad%g2rt(kt+i)=drdr+2.0*rdr/rv1+&
170 7729260 : (rdtt+rdt/tant1+rdff/sint2)/rv2
171 :
172 7729260 : dzdr = ((rdru-rdrd)*ro- (rou-rod)*rdr)/ro2
173 :
174 : !-- > dzdtr,dzdfs vanish by definition.
175 :
176 7729260 : dzdtr = 0.0
177 7729260 : dzdfs = 0.0
178 :
179 7729260 : grad%gggrt(kt+i) = grr*dagrr + grt*dagrt + grf*dagrf
180 :
181 7729260 : grad%gzgr(kt+i) = dzdr*grr + dzdtr*grt + dzdfs*grf
182 :
183 7729260 : grru = rdru
184 7729260 : grtu = rdtu/rv1
185 7729260 : grfu = rdfu/rvsin1
186 :
187 7729260 : grad%agru(kt+i) = sqrt(grru**2+grtu**2+grfu**2)
188 :
189 : dagrru = (rdru*rdrru*rv3+rdtu* (rdrtu*rv1-rdtu)+&
190 7729260 : rdfu* (rdrfu*rv1-rdfu)/sint2)/grad%agru(kt+i)/rv3
191 :
192 : dagrtu = (rdru*rdrtu*rv2+rdtu*rdttu+&
193 7729260 : rdfu* (-rdfu/tant1+rdtfu)/sint2)/ (grad%agru(kt+i)*rv3)
194 :
195 : dagrfu = (rdru*rdrfu*rv2+rdtu*rdtfu+rdfu*rdffu/sint2)/&
196 7729260 : (grad%agru(kt+i)*rv3*sint1)
197 :
198 : grad%g2ru(kt+i) = rdrru + 2.e0*rdru/rv1 +&
199 7729260 : (rdttu+rdtu/tant1+rdffu/sint2)/rv2
200 :
201 7729260 : grad%gggru(kt+i) = grru*dagrru + grtu*dagrtu + grfu*dagrfu
202 :
203 7729260 : grad%grgru(kt+i) = grr*grru + grt*grtu + grf*grfu
204 :
205 7729260 : grrd = rdrd
206 7729260 : grtd = rdtd/rv1
207 7729260 : grfd = rdfd/rvsin1
208 :
209 7729260 : grad%agrd(kt+i) = sqrt(grrd**2+grtd**2+grfd**2)
210 :
211 : dagrrd = (rdrd*rdrrd*rv3+rdtd* (rdrtd*rv1-rdtd)+&
212 7729260 : rdfd* (rdrfd*rv1-rdfd)/sint2)/grad%agrd(kt+i)/rv3
213 :
214 : dagrtd = (rdrd*rdrtd*rv2+rdtd*rdttd+&
215 7729260 : rdfd* (-rdfd/tant1+rdtfd)/sint2)/ (grad%agrd(kt+i)*rv3)
216 :
217 : dagrfd = (rdrd*rdrfd*rv2+rdtd*rdtfd+rdfd*rdffd/sint2)/&
218 7729260 : (grad%agrd(kt+i)*rv3*sint1)
219 :
220 : grad%g2rd(kt+i) = rdrrd + 2*rdrd/rv1 +&
221 7729260 : (rdttd+rdtd/tant1+rdffd/sint2)/rv2
222 :
223 7729260 : grad%gggrd(kt+i) = grrd*dagrrd + grtd*dagrtd + grfd*dagrfd
224 :
225 7766442 : grad%grgrd(kt+i) = grr*grrd + grt*grtd + grf*grfd
226 :
227 : ENDDO points_1
228 :
229 : ELSE
230 :
231 21449427 : points_2 : DO i = 1,nsp
232 :
233 21324888 : grad%agrt(kt+i) = 0.0
234 21324888 : grad%agru(kt+i) = 0.0
235 21324888 : grad%agrd(kt+i) = 0.0
236 21324888 : grad%g2rt(kt+i) = 0.0
237 21324888 : grad%g2ru(kt+i) = 0.0
238 21324888 : grad%g2rd(kt+i) = 0.0
239 21324888 : grad%gggrt(kt+i) = 0.0
240 21324888 : grad%gggru(kt+i) = 0.0
241 21324888 : grad%gggrd(kt+i) = 0.0
242 21324888 : grad%grgru(kt+i) = 0.0
243 21324888 : grad%grgrd(kt+i) = 0.0
244 21324888 : grad%gzgr(kt+i) = 0.0
245 :
246 21324888 : ro = rh(i,1) + rh(i,jspins)
247 :
248 21324888 : IF (ro<chsml) CYCLE points_2
249 :
250 21324888 : sint1 = sin(thet(i))
251 21324888 : sint2 = sint1**2
252 21324888 : tant1 = tan(thet(i))
253 21324888 : rvsin1 = rv1*sint1
254 :
255 21324888 : rou = rh(i,1)
256 21324888 : rdru = rhdr(i,1)
257 21324888 : rdtu = rhdt(i,1)
258 21324888 : rdfu = rhdf(i,1)
259 21324888 : rdrru = rhdrr(i,1)
260 21324888 : rdttu = rhdtt(i,1)
261 21324888 : rdffu = rhdff(i,1)
262 21324888 : rdtfu = rhdtf(i,1)
263 21324888 : rdrtu = rhdrt(i,1)
264 21324888 : rdrfu = rhdrf(i,1)
265 :
266 21324888 : rod = rh(i,jspins)
267 21324888 : rdrd = rhdr(i,jspins)
268 21324888 : rdtd = rhdt(i,jspins)
269 21324888 : rdfd = rhdf(i,jspins)
270 21324888 : rdrrd = rhdrr(i,jspins)
271 21324888 : rdttd = rhdtt(i,jspins)
272 21324888 : rdffd = rhdff(i,jspins)
273 21324888 : rdtfd = rhdtf(i,jspins)
274 21324888 : rdrtd = rhdrt(i,jspins)
275 21324888 : rdrfd = rhdrf(i,jspins)
276 :
277 21324888 : rdr = rdru + rdrd
278 21324888 : rdt = rdtu + rdtd
279 21324888 : rdf = rdfu + rdfd
280 21324888 : drdr = rdrru + rdrrd
281 21324888 : rdtt = rdttu + rdttd
282 21324888 : rdff = rdffu + rdffd
283 21324888 : rdtf = rdtfu + rdtfd
284 21324888 : rdrt = rdrtu + rdrtd
285 21324888 : rdrf = rdrfu + rdrfd
286 :
287 21324888 : ro2 = ro**2
288 :
289 21324888 : grr = rdr
290 21324888 : grt = rdt/rv1
291 21324888 : grf = rdf/rvsin1
292 :
293 21324888 : grad%agrt(kt+i) = sqrt(grr**2+grt**2+grf**2)
294 :
295 21324888 : IF (grad%agrt(kt+i)<chsml) CYCLE points_2
296 :
297 21324888 : dagrr = (rdr*drdr*rv3+rdt* (rdrt*rv1-rdt)+ rdf* (rdrf*rv1-rdf)/sint2)/grad%agrt(kt+i)/rv3
298 :
299 21324888 : dagrt =(rdr*rdrt*rv2+rdt*rdtt+rdf* (-rdf/tant1+rdtf)/sint2)/ (grad%agrt(kt+i)*rv3)
300 :
301 21324888 : dagrf = (rdr*rdrf*rv2+rdt*rdtf+rdf*rdff/sint2)/(grad%agrt(kt+i)*rv3*sint1)
302 :
303 21324888 : grad%g2rt(kt+i)= drdr+2.0*rdr/rv1+(rdtt+rdt/tant1+rdff/sint2)/rv2
304 :
305 21324888 : dzdr = ((rdru-rdrd)*ro- (rou-rod)*rdr)/ro2
306 :
307 : ! dzdtr,dzdfs vanish by definition.
308 21324888 : dzdtr = 0.0
309 21324888 : dzdfs = 0.0
310 :
311 21324888 : grad%gggrt(kt+i) = grr*dagrr + grt*dagrt + grf*dagrf
312 :
313 21324888 : grad%gzgr(kt+i) = dzdr*grr + dzdtr*grt + dzdfs*grf
314 :
315 21324888 : grru = rdru
316 21324888 : grtu = rdtu/rv1
317 21324888 : grfu = rdfu/rvsin1
318 :
319 21324888 : grad%agru(kt+i) = sqrt(grru**2+grtu**2+grfu**2)
320 :
321 : dagrru = (rdru*rdrru*rv3+rdtu* (rdrtu*rv1-rdtu)+&
322 21324888 : rdfu* (rdrfu*rv1-rdfu)/sint2)/grad%agru(kt+i)/rv3
323 :
324 21324888 : dagrtu = (rdru*rdrtu*rv2+rdtu*rdttu+ rdfu* (-rdfu/tant1+rdtfu)/sint2)/ (grad%agru(kt+i)*rv3)
325 :
326 21324888 : dagrfu = (rdru*rdrfu*rv2+rdtu*rdtfu+rdfu*rdffu/sint2)/ (grad%agru(kt+i)*rv3*sint1)
327 :
328 21324888 : grad%g2ru(kt+i) = rdrru + 2.e0*rdru/rv1 + (rdttu+rdtu/tant1+rdffu/sint2)/rv2
329 :
330 21324888 : grad%gggru(kt+i) = grru*dagrru + grtu*dagrtu + grfu*dagrfu
331 :
332 21324888 : grad%grgru(kt+i) = grr*grru + grt*grtu + grf*grfu
333 :
334 21324888 : grrd = rdrd
335 21324888 : grtd = rdtd/rv1
336 21324888 : grfd = rdfd/rvsin1
337 :
338 21324888 : grad%agrd(kt+i) = sqrt(grrd**2+grtd**2+grfd**2)
339 :
340 21324888 : dagrrd = (rdrd*rdrrd*rv3+rdtd* (rdrtd*rv1-rdtd)+ rdfd* (rdrfd*rv1-rdfd)/sint2)/grad%agrd(kt+i)/rv3
341 :
342 21324888 : dagrtd = (rdrd*rdrtd*rv2+rdtd*rdttd+ rdfd* (-rdfd/tant1+rdtfd)/sint2)/ (grad%agrd(kt+i)*rv3)
343 :
344 21324888 : dagrfd = (rdrd*rdrfd*rv2+rdtd*rdtfd+rdfd*rdffd/sint2)/ (grad%agrd(kt+i)*rv3*sint1)
345 :
346 21324888 : grad%g2rd(kt+i) = rdrrd + 2*rdrd/rv1 + (rdttd+rdtd/tant1+rdffd/sint2)/rv2
347 :
348 21324888 : grad%gggrd(kt+i) = grrd*dagrrd + grtd*dagrtd + grfd*dagrfd
349 :
350 21449427 : grad%grgrd(kt+i) = grr*grrd + grt*grtd + grf*grfd
351 :
352 : ENDDO points_2
353 :
354 : ENDIF
355 : ENDIF
356 : END SUBROUTINE mkgylm
357 : END MODULE m_mkgylm
|
