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 : !-----------------------------------------------------------------------------
8 : ! Using the components and derivatives of a charge density rho in spherical
9 : ! coordinates on the real space grid, make the following quantaties:
10 : !
11 : ! gr(js): [partial_r (rho),partial_theta (rho),partial_phi (rho)]
12 : ! sigma: |grad(rho)|^2 for jspins==1, otherwise three components, namely
13 : ! |grad(rho_up)|^2,grad(rho_up)*grad(rho_down),|grad(rho_down)|^2
14 : ! laplace(js): laplace(rho_js)
15 : !
16 : ! and these older components:
17 : ! agrt/u/d: |grad(rho)| for total density/spin-up/spin-down
18 : ! g2rt/u/d: laplace(rho)
19 : ! gggrt/u/d: (grad(rho))*(grad(|grad(rho)|)) [scalar product]
20 : ! gzgr: (grad(zeta))*(grad(rho)) for zeta=(rho_up-rho_down)/rho
21 : !
22 : ! which are used to calculate gradient contributions to the xc potential and
23 : ! energy.
24 : !
25 : ! rh is rho, rhd[i][j] are the partial derivatives along one/two directions.
26 : ! rv and thet are the radius and polar angles and t/f in derivatives stand for
27 : ! theta/phi.
28 : !
29 : ! Modified so only allocated old quantities require calculations. A.N. 2019
30 : !
31 : ! Quantities fo libxc are calculated as well. D.W./M.R. 2018
32 : !
33 : ! Original script by T.A. 1996
34 : !-----------------------------------------------------------------------------
35 : CONTAINS
36 338811 : SUBROUTINE mkgylm(jspins,rv,thet,nsp, rh,rhdr,rhdt,rhdf,rhdrr,rhdtt,rhdff, rhdtf,rhdrt,rhdrf,grad,kt)
37 : USE m_types
38 : IMPLICIT NONE
39 : REAL, INTENT (IN) :: rv
40 : REAL, INTENT (IN) :: thet(:)
41 : REAL, INTENT (IN) :: rh(:,:)
42 : REAL, INTENT (IN) :: rhdr(:,:),rhdt(:,:),rhdf(:,:)
43 : REAL, INTENT (IN) :: rhdrr(:,:),rhdtt(:,:),rhdff(:,:)
44 : REAL, INTENT (IN) :: rhdtf(:,:),rhdrf(:,:),rhdrt(:,:)
45 : TYPE(t_gradients), INTENT(INOUT) :: grad
46 : INTEGER, INTENT(IN) :: jspins,nsp
47 : INTEGER, INTENT(IN) :: kt !index of first point to use in gradients
48 :
49 : REAL dagrf,dagrfd,dagrfu,dagrr,dagrrd,dagrru,dagrt, &
50 : dagrtd,dagrtu,drdr,dzdfs,dzdr,dzdtr,grf,grfd,grfu,grr, &
51 : grrd,grru,grt,grtd,grtu,rdf,rdfd,rdff,rdffd,rdffu,rdfu, &
52 : rdr,rdrd,rdrf,rdrfd,rdrfu,rdrrd,rdrru,rdrt,rdrtd,rdrtu, &
53 : rdru,rdt,rdtd,rdtf,rdtfd,rdtfu,rdtt,rdttd,rdttu,rdtu, &
54 : ro,ro2,rod,rou,rv1,rv2,rv3,rvsin1,sint1,sint2,tant1
55 : REAL, PARAMETER :: chsml = 1.e-10
56 : INTEGER i,js,fac
57 :
58 338811 : rv1 = rv
59 338811 : rv2 = rv1**2
60 338811 : rv3 = rv1**3
61 :
62 : ! Gradients for libxc calculations.
63 338811 : IF (ALLOCATED(grad%gr)) THEN
64 42252 : DO js=1,jspins
65 5521002 : DO i=1,nsp
66 21944601 : grad%gr(:,kt+i,js)=[rhdr(i,js),rhdt(i,js)/rv1,rhdf(i,js)/(rv1*sin(thet(i)))]
67 : END DO
68 : END DO
69 : ! Use contracted gradients only for libxc.
70 12651 : IF (ALLOCATED(grad%sigma)) THEN
71 4217 : IF (jspins==1) THEN
72 290277 : DO i=1,nsp
73 1158627 : grad%sigma(1,kt+i)= dot_PRODUCT(grad%gr(:,kt+i,1),grad%gr(:,kt+i,1))
74 : END DO
75 : ELSE
76 579690 : DO i=1,nsp
77 2305200 : grad%sigma(1,kt+i)= dot_PRODUCT(grad%gr(:,kt+i,1),grad%gr(:,kt+i,1))
78 2305200 : grad%sigma(2,kt+i)= dot_PRODUCT(grad%gr(:,kt+i,1),grad%gr(:,kt+i,2))
79 2308590 : grad%sigma(3,kt+i)= dot_PRODUCT(grad%gr(:,kt+i,2),grad%gr(:,kt+i,2))
80 : END DO
81 : END IF
82 : END IF
83 12651 : IF (ALLOCATED(grad%laplace)) THEN
84 : !Lapacian of density
85 : !fac=MERGE(2,1,jspins==1)
86 : fac=1
87 11824 : DO js=1,jspins
88 1453874 : DO i=1,nsp
89 : grad%laplace(kt+i,js) = (rhdrr(i,js) + 2.e0*rhdr(i,js)/rv1 +&
90 1449657 : (rhdtt(i,js)+rhdt(i,js)/TAN(thet(i))+rhdff(i,js)/SIN(thet(i))**2)/rv2)/fac
91 : ENDDO
92 : ENDDO
93 : ENDIF
94 : ! Cartesian components of the gradient for sourcefree calculations. TODO: Need phi here as well.
95 : !IF (ALLOCATED(grad%grxyz)) THEN
96 : ! grad%grxyz=SIN(th)*COS(ph)*grad%gr(1,kt+k,1) + COS(th)*COS(ph)*grad%gr(2,kt+k,1)/r - SIN(ph)*grad%gr(3,kt+k,1)/(r*SIN(th))
97 : !END IF
98 12651 : RETURN ! Do not calculate arrays for in-build GGA.
99 : END IF
100 :
101 : IF (ALLOCATED(grad%gr)) THEN
102 :
103 : END IF
104 : ! Old code for in-build xcpots
105 326160 : IF(allocated(grad%agrt)) THEN
106 326160 : IF (jspins==1) THEN
107 :
108 34433706 : points_1 : DO i = 1,nsp
109 :
110 34262388 : grad%agrt(kt+i) = 0.0
111 34262388 : grad%agru(kt+i) = 0.0
112 34262388 : grad%agrd(kt+i) = 0.0
113 34262388 : grad%g2rt(kt+i) = 0.0
114 34262388 : grad%g2ru(kt+i) = 0.0
115 34262388 : grad%g2rd(kt+i) = 0.0
116 34262388 : grad%gggrt(kt+i) = 0.0
117 34262388 : grad%gggru(kt+i) = 0.0
118 34262388 : grad%gggrd(kt+i) = 0.0
119 34262388 : grad%grgru(kt+i) = 0.0
120 34262388 : grad%grgrd(kt+i) = 0.0
121 34262388 : grad%gzgr(kt+i) = 0.0
122 :
123 34262388 : ro = rh(i,1)
124 :
125 34262388 : IF (ro<chsml) CYCLE points_1
126 34262388 : sint1 = sin(thet(i))
127 34262388 : sint2 = sint1**2
128 34262388 : tant1 = tan(thet(i))
129 34262388 : rvsin1 = rv1*sint1
130 :
131 34262388 : rou = ro/2
132 34262388 : rdru = rhdr(i,1)/2
133 34262388 : rdtu = rhdt(i,1)/2
134 34262388 : rdfu = rhdf(i,1)/2
135 34262388 : rdrru = rhdrr(i,1)/2
136 34262388 : rdttu = rhdtt(i,1)/2
137 34262388 : rdffu = rhdff(i,1)/2
138 34262388 : rdtfu = rhdtf(i,1)/2
139 34262388 : rdrtu = rhdrt(i,1)/2
140 34262388 : rdrfu = rhdrf(i,1)/2
141 :
142 34262388 : rod = rou
143 34262388 : rdrd = rdru
144 34262388 : rdtd = rdtu
145 34262388 : rdfd = rdfu
146 34262388 : rdrrd = rdrru
147 34262388 : rdttd = rdttu
148 34262388 : rdffd = rdffu
149 34262388 : rdtfd = rdtfu
150 34262388 : rdrtd = rdrtu
151 34262388 : rdrfd = rdrfu
152 :
153 34262388 : rdr = rdru + rdrd
154 34262388 : rdt = rdtu + rdtd
155 34262388 : rdf = rdfu + rdfd
156 34262388 : drdr = rdrru + rdrrd
157 34262388 : rdtt = rdttu + rdttd
158 34262388 : rdff = rdffu + rdffd
159 34262388 : rdtf = rdtfu + rdtfd
160 34262388 : rdrt = rdrtu + rdrtd
161 34262388 : rdrf = rdrfu + rdrfd
162 :
163 34262388 : ro2 = ro**2
164 :
165 34262388 : grr = rdr
166 34262388 : grt = rdt/rv1
167 34262388 : grf = rdf/rvsin1
168 :
169 34262388 : grad%agrt(kt+i) = sqrt(grr**2+grt**2+grf**2)
170 :
171 34262388 : IF (grad%agrt(kt+i)<chsml) CYCLE points_1
172 :
173 : dagrr = (rdr*drdr*rv3+rdt* (rdrt*rv1-rdt)+ &
174 34262388 : rdf* (rdrf*rv1-rdf)/sint2)/grad%agrt(kt+i)/rv3
175 :
176 : dagrt =(rdr*rdrt*rv2+rdt*rdtt+rdf* (-rdf/tant1+rdtf)/sint2)/&
177 34262388 : (grad%agrt(kt+i)*rv3)
178 :
179 : dagrf = (rdr*rdrf*rv2+rdt*rdtf+rdf*rdff/sint2)/&
180 34262388 : (grad%agrt(kt+i)*rv3*sint1)
181 :
182 : grad%g2rt(kt+i)=drdr+2.0*rdr/rv1+&
183 34262388 : (rdtt+rdt/tant1+rdff/sint2)/rv2
184 :
185 34262388 : dzdr = ((rdru-rdrd)*ro- (rou-rod)*rdr)/ro2
186 :
187 : !-- > dzdtr,dzdfs vanish by definition.
188 :
189 34262388 : dzdtr = 0.0
190 34262388 : dzdfs = 0.0
191 :
192 34262388 : grad%gggrt(kt+i) = grr*dagrr + grt*dagrt + grf*dagrf
193 :
194 34262388 : grad%gzgr(kt+i) = dzdr*grr + dzdtr*grt + dzdfs*grf
195 :
196 34262388 : grru = rdru
197 34262388 : grtu = rdtu/rv1
198 34262388 : grfu = rdfu/rvsin1
199 :
200 34262388 : grad%agru(kt+i) = sqrt(grru**2+grtu**2+grfu**2)
201 :
202 : dagrru = (rdru*rdrru*rv3+rdtu* (rdrtu*rv1-rdtu)+&
203 34262388 : rdfu* (rdrfu*rv1-rdfu)/sint2)/grad%agru(kt+i)/rv3
204 :
205 : dagrtu = (rdru*rdrtu*rv2+rdtu*rdttu+&
206 34262388 : rdfu* (-rdfu/tant1+rdtfu)/sint2)/ (grad%agru(kt+i)*rv3)
207 :
208 : dagrfu = (rdru*rdrfu*rv2+rdtu*rdtfu+rdfu*rdffu/sint2)/&
209 34262388 : (grad%agru(kt+i)*rv3*sint1)
210 :
211 : grad%g2ru(kt+i) = rdrru + 2.e0*rdru/rv1 +&
212 34262388 : (rdttu+rdtu/tant1+rdffu/sint2)/rv2
213 :
214 34262388 : grad%gggru(kt+i) = grru*dagrru + grtu*dagrtu + grfu*dagrfu
215 :
216 34262388 : grad%grgru(kt+i) = grr*grru + grt*grtu + grf*grfu
217 :
218 34262388 : grrd = rdrd
219 34262388 : grtd = rdtd/rv1
220 34262388 : grfd = rdfd/rvsin1
221 :
222 34262388 : grad%agrd(kt+i) = sqrt(grrd**2+grtd**2+grfd**2)
223 :
224 : dagrrd = (rdrd*rdrrd*rv3+rdtd* (rdrtd*rv1-rdtd)+&
225 34262388 : rdfd* (rdrfd*rv1-rdfd)/sint2)/grad%agrd(kt+i)/rv3
226 :
227 : dagrtd = (rdrd*rdrtd*rv2+rdtd*rdttd+&
228 34262388 : rdfd* (-rdfd/tant1+rdtfd)/sint2)/ (grad%agrd(kt+i)*rv3)
229 :
230 : dagrfd = (rdrd*rdrfd*rv2+rdtd*rdtfd+rdfd*rdffd/sint2)/&
231 34262388 : (grad%agrd(kt+i)*rv3*sint1)
232 :
233 : grad%g2rd(kt+i) = rdrrd + 2*rdrd/rv1 +&
234 34262388 : (rdttd+rdtd/tant1+rdffd/sint2)/rv2
235 :
236 34262388 : grad%gggrd(kt+i) = grrd*dagrrd + grtd*dagrtd + grfd*dagrfd
237 :
238 34433706 : grad%grgrd(kt+i) = grr*grrd + grt*grtd + grf*grfd
239 :
240 : ENDDO points_1
241 :
242 : ELSE
243 :
244 29781844 : points_2 : DO i = 1,nsp
245 :
246 29627002 : grad%agrt(kt+i) = 0.0
247 29627002 : grad%agru(kt+i) = 0.0
248 29627002 : grad%agrd(kt+i) = 0.0
249 29627002 : grad%g2rt(kt+i) = 0.0
250 29627002 : grad%g2ru(kt+i) = 0.0
251 29627002 : grad%g2rd(kt+i) = 0.0
252 29627002 : grad%gggrt(kt+i) = 0.0
253 29627002 : grad%gggru(kt+i) = 0.0
254 29627002 : grad%gggrd(kt+i) = 0.0
255 29627002 : grad%grgru(kt+i) = 0.0
256 29627002 : grad%grgrd(kt+i) = 0.0
257 29627002 : grad%gzgr(kt+i) = 0.0
258 :
259 29627002 : ro = rh(i,1) + rh(i,jspins)
260 :
261 29627002 : IF (ro<chsml) CYCLE points_2
262 :
263 29627002 : sint1 = sin(thet(i))
264 29627002 : sint2 = sint1**2
265 29627002 : tant1 = tan(thet(i))
266 29627002 : rvsin1 = rv1*sint1
267 :
268 29627002 : rou = rh(i,1)
269 29627002 : rdru = rhdr(i,1)
270 29627002 : rdtu = rhdt(i,1)
271 29627002 : rdfu = rhdf(i,1)
272 29627002 : rdrru = rhdrr(i,1)
273 29627002 : rdttu = rhdtt(i,1)
274 29627002 : rdffu = rhdff(i,1)
275 29627002 : rdtfu = rhdtf(i,1)
276 29627002 : rdrtu = rhdrt(i,1)
277 29627002 : rdrfu = rhdrf(i,1)
278 :
279 29627002 : rod = rh(i,jspins)
280 29627002 : rdrd = rhdr(i,jspins)
281 29627002 : rdtd = rhdt(i,jspins)
282 29627002 : rdfd = rhdf(i,jspins)
283 29627002 : rdrrd = rhdrr(i,jspins)
284 29627002 : rdttd = rhdtt(i,jspins)
285 29627002 : rdffd = rhdff(i,jspins)
286 29627002 : rdtfd = rhdtf(i,jspins)
287 29627002 : rdrtd = rhdrt(i,jspins)
288 29627002 : rdrfd = rhdrf(i,jspins)
289 :
290 29627002 : rdr = rdru + rdrd
291 29627002 : rdt = rdtu + rdtd
292 29627002 : rdf = rdfu + rdfd
293 29627002 : drdr = rdrru + rdrrd
294 29627002 : rdtt = rdttu + rdttd
295 29627002 : rdff = rdffu + rdffd
296 29627002 : rdtf = rdtfu + rdtfd
297 29627002 : rdrt = rdrtu + rdrtd
298 29627002 : rdrf = rdrfu + rdrfd
299 :
300 29627002 : ro2 = ro**2
301 :
302 29627002 : grr = rdr
303 29627002 : grt = rdt/rv1
304 29627002 : grf = rdf/rvsin1
305 :
306 29627002 : grad%agrt(kt+i) = sqrt(grr**2+grt**2+grf**2)
307 :
308 29627002 : IF (grad%agrt(kt+i)<chsml) CYCLE points_2
309 :
310 29627002 : dagrr = (rdr*drdr*rv3+rdt* (rdrt*rv1-rdt)+ rdf* (rdrf*rv1-rdf)/sint2)/grad%agrt(kt+i)/rv3
311 :
312 29627002 : dagrt =(rdr*rdrt*rv2+rdt*rdtt+rdf* (-rdf/tant1+rdtf)/sint2)/ (grad%agrt(kt+i)*rv3)
313 :
314 29627002 : dagrf = (rdr*rdrf*rv2+rdt*rdtf+rdf*rdff/sint2)/(grad%agrt(kt+i)*rv3*sint1)
315 :
316 29627002 : grad%g2rt(kt+i)= drdr+2.0*rdr/rv1+(rdtt+rdt/tant1+rdff/sint2)/rv2
317 :
318 29627002 : dzdr = ((rdru-rdrd)*ro- (rou-rod)*rdr)/ro2
319 :
320 : ! dzdtr,dzdfs vanish by definition.
321 29627002 : dzdtr = 0.0
322 29627002 : dzdfs = 0.0
323 :
324 29627002 : grad%gggrt(kt+i) = grr*dagrr + grt*dagrt + grf*dagrf
325 :
326 29627002 : grad%gzgr(kt+i) = dzdr*grr + dzdtr*grt + dzdfs*grf
327 :
328 29627002 : grru = rdru
329 29627002 : grtu = rdtu/rv1
330 29627002 : grfu = rdfu/rvsin1
331 :
332 29627002 : grad%agru(kt+i) = sqrt(grru**2+grtu**2+grfu**2)
333 :
334 : dagrru = (rdru*rdrru*rv3+rdtu* (rdrtu*rv1-rdtu)+&
335 29627002 : rdfu* (rdrfu*rv1-rdfu)/sint2)/grad%agru(kt+i)/rv3
336 :
337 29627002 : dagrtu = (rdru*rdrtu*rv2+rdtu*rdttu+ rdfu* (-rdfu/tant1+rdtfu)/sint2)/ (grad%agru(kt+i)*rv3)
338 :
339 29627002 : dagrfu = (rdru*rdrfu*rv2+rdtu*rdtfu+rdfu*rdffu/sint2)/ (grad%agru(kt+i)*rv3*sint1)
340 :
341 29627002 : grad%g2ru(kt+i) = rdrru + 2.e0*rdru/rv1 + (rdttu+rdtu/tant1+rdffu/sint2)/rv2
342 :
343 29627002 : grad%gggru(kt+i) = grru*dagrru + grtu*dagrtu + grfu*dagrfu
344 :
345 29627002 : grad%grgru(kt+i) = grr*grru + grt*grtu + grf*grfu
346 :
347 29627002 : grrd = rdrd
348 29627002 : grtd = rdtd/rv1
349 29627002 : grfd = rdfd/rvsin1
350 :
351 29627002 : grad%agrd(kt+i) = sqrt(grrd**2+grtd**2+grfd**2)
352 :
353 29627002 : dagrrd = (rdrd*rdrrd*rv3+rdtd* (rdrtd*rv1-rdtd)+ rdfd* (rdrfd*rv1-rdfd)/sint2)/grad%agrd(kt+i)/rv3
354 :
355 29627002 : dagrtd = (rdrd*rdrtd*rv2+rdtd*rdttd+ rdfd* (-rdfd/tant1+rdtfd)/sint2)/ (grad%agrd(kt+i)*rv3)
356 :
357 29627002 : dagrfd = (rdrd*rdrfd*rv2+rdtd*rdtfd+rdfd*rdffd/sint2)/ (grad%agrd(kt+i)*rv3*sint1)
358 :
359 29627002 : grad%g2rd(kt+i) = rdrrd + 2*rdrd/rv1 + (rdttd+rdtd/tant1+rdffd/sint2)/rv2
360 :
361 29627002 : grad%gggrd(kt+i) = grrd*dagrrd + grtd*dagrtd + grfd*dagrfd
362 :
363 29781844 : grad%grgrd(kt+i) = grr*grrd + grt*grtd + grf*grfd
364 :
365 : ENDDO points_2
366 :
367 : ENDIF
368 : ENDIF
369 : END SUBROUTINE mkgylm
370 : END MODULE m_mkgylm
|
