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_mt_tofrom_grid
7 : USE m_types
8 : PRIVATE
9 : REAL, PARAMETER :: d_15 = 1.e-15
10 : REAL, ALLOCATABLE :: ylh(:, :, :), ylht(:, :, :), ylhtt(:, :, :)
11 : REAL, ALLOCATABLE :: ylhf(:, :, :), ylhff(:, :, :), ylhtf(:, :, :)
12 : REAL, ALLOCATABLE :: wt(:), rx(:, :), thet(:), phi(:)
13 : PUBLIC :: init_mt_grid, mt_to_grid, mt_from_grid, finish_mt_grid
14 : CONTAINS
15 732 : SUBROUTINE init_mt_grid(jspins, atoms, sphhar, dograds, sym, thout, phout)
16 : USE m_gaussp
17 : USE m_lhglptg
18 : USE m_lhglpts
19 : IMPLICIT NONE
20 : INTEGER, INTENT(IN) :: jspins
21 : TYPE(t_atoms), INTENT(IN) :: atoms
22 : TYPE(t_sphhar), INTENT(IN) :: sphhar
23 : LOGICAL, INTENT(IN) :: dograds
24 : TYPE(t_sym), INTENT(IN) :: sym
25 : REAL, INTENT(OUT), OPTIONAL :: thout(:)
26 : REAL, INTENT(OUT), OPTIONAL :: phout(:)
27 :
28 732 : call timestart("init_mt_grid")
29 : ! generate nspd points on a sherical shell with radius 1.0
30 : ! angular mesh equidistant in phi,
31 : ! theta are zeros of the legendre polynomials
32 6588 : ALLOCATE (wt(atoms%nsp()), rx(3, atoms%nsp()), thet(atoms%nsp()), phi(atoms%nsp()))
33 732 : CALL gaussp(atoms%lmaxd, rx, wt)
34 : ! generate the lattice harmonics on the angular mesh
35 3660 : ALLOCATE (ylh(atoms%nsp(), 0:sphhar%nlhd, sphhar%ntypsd))
36 732 : IF (dograds) THEN
37 2830 : ALLOCATE (ylht, MOLD=ylh)
38 2264 : ALLOCATE (ylhtt, MOLD=ylh)
39 2264 : ALLOCATE (ylhf, MOLD=ylh)
40 2264 : ALLOCATE (ylhff, MOLD=ylh)
41 2264 : ALLOCATE (ylhtf, MOLD=ylh)
42 :
43 : CALL lhglptg(sphhar, atoms, rx, atoms%nsp(), dograds, sym, &
44 566 : ylh, thet, phi, ylht, ylhtt, ylhf, ylhff, ylhtf)
45 566 : IF (PRESENT(thout)) THEN
46 0 : thout=thet
47 0 : phout=phi
48 : END IF
49 :
50 : ELSE
51 166 : CALL lhglpts(sphhar, atoms, rx, atoms%nsp(), sym, ylh)
52 : END IF
53 732 : call timestop("init_mt_grid")
54 732 : END SUBROUTINE init_mt_grid
55 :
56 692 : SUBROUTINE mt_to_grid(dograds, jspins, atoms, sym,sphhar,rotch, den_mt, n, noco ,grad, ch)
57 : USE m_grdchlh
58 : USE m_mkgylm
59 : IMPLICIT NONE
60 : LOGICAL, INTENT(IN) :: dograds
61 : TYPE(t_atoms), INTENT(IN) :: atoms
62 : TYPE(t_sym), INTENT(IN) :: sym
63 : LOGICAL, INTENT(IN) :: rotch
64 : TYPE(t_sphhar), INTENT(IN) :: sphhar
65 : REAL, INTENT(IN) :: den_mt(:, 0:, :)
66 : INTEGER, INTENT(IN) :: n, jspins
67 : REAL, INTENT(OUT), OPTIONAL :: ch(:, :)
68 : TYPE(t_gradients), INTENT(INOUT):: grad
69 : TYPE(t_noco), INTENT(IN) :: noco
70 : REAL :: dentot
71 :
72 :
73 :
74 : REAL :: rho_11,rho_22,rho_21r,rho_21i,mx,my,mz,magmom
75 : REAL :: rhotot,rho_up,rho_down
76 692 : REAL, ALLOCATABLE :: chlh(:, :, :), chlhdr(:, :, :), chlhdrr(:, :, :)
77 692 : REAL, ALLOCATABLE :: chdr(:, :), chdt(:, :), chdf(:, :), ch_tmp(:, :),ch_calc(:,:)
78 692 : REAL, ALLOCATABLE :: chdrr(:, :), chdtt(:, :), chdff(:, :), chdtf(:, :)
79 692 : REAL, ALLOCATABLE :: chdrt(:, :), chdrf(:, :)
80 692 : REAL, ALLOCATABLE :: drm(:,:), drrm(:,:), mm(:,:)
81 692 : REAL, ALLOCATABLE :: chlhtot(:),chlhdrtot(:),chlhdrrtot(:)
82 : INTEGER:: nd, lh, js, jr, kt, k, nsp,j,i,jspV
83 :
84 692 : call timestart("mt_to_grid")
85 : !This snippet is crucial to determine over which spins (Only diagonals in colinear case or also off diags in non colin case.)
86 1987 : IF (any(noco%l_unrestrictMT)) THEN
87 : jspV=4
88 : ELSE
89 607 : jspV=jspins
90 : END IF
91 :
92 692 : nd = sym%ntypsy(atoms%firstAtom(n))
93 692 : nsp = atoms%nsp()
94 :
95 : !General Allocations
96 3460 : ALLOCATE (chlh(atoms%jmtd, 0:sphhar%nlhd, jspV))
97 4844 : ALLOCATE (ch_tmp(nsp, jspV),ch_calc(nsp*atoms%jmtd, jspV))
98 :
99 : !Allocations in dograds case
100 692 : IF (dograds) THEN
101 0 : ALLOCATE (chdr(nsp, jspV), chdt(nsp, jspV), chdf(nsp, jspV), chdrr(nsp, jspV), &
102 0 : chdtt(nsp, jspV), chdff(nsp, jspV), chdtf(nsp, jspV), chdrt(nsp, jspV), &
103 9785 : chdrf(nsp, jspV))
104 2060 : ALLOCATE (chlhdr(atoms%jmtd, 0:sphhar%nlhd, jspV))
105 2060 : ALLOCATE (chlhdrr(atoms%jmtd, 0:sphhar%nlhd, jspV))
106 : ENDIF
107 :
108 : !Allocations in mtNoco case
109 1987 : IF (any(noco%l_unrestrictMT)) THEN
110 : !General Noco Allocations
111 340 : ALLOCATE(mm(atoms%jmtd, 0:sphhar%nlhd))
112 :
113 : !Allocations in case one uses e.g. GGA with mtNoco
114 85 : IF (dograds) THEN
115 0 : ALLOCATE(drm(atoms%jmtd,0:sphhar%nlhd),drrm(atoms%jmtd, 0:sphhar%nlhd))
116 0 : ALLOCATE(chlhtot(0:sphhar%nlhd),chlhdrtot(0:sphhar%nlhd),chlhdrrtot(0:sphhar%nlhd))
117 : END IF
118 : END IF
119 :
120 : !Calc magnetization (This is necessary only in mtNoco case)
121 5160667 : IF(any(noco%l_unrestrictMT)) mm(:,:)=SQRT((0.5*(den_mt(:,:,1)-den_mt(:,:,2)))**2+4*den_mt(:,:,3)**2+4*den_mt(:,:,4)**2)
122 :
123 :
124 : !Loop to calculate chlh and necessary gradients (if needed)
125 21724 : DO lh = 0, sphhar%nlh(nd)
126 : ! calculates gradients of radial charge densities of l=> 0.
127 : ! rho*ylh/r**2 is charge density. chlh=rho/r**2.
128 : ! charge density=sum(chlh*ylh).
129 : ! chlhdr=d(chlh)/dr, chlhdrr=dd(chlh)/drr.
130 :
131 : !Scaling of the magnetic moments in the same way the charge density is scaled in chlh.
132 5197931 : IF(any(noco%l_unrestrictMT)) mm(:,lh)=0.5*mm(:,lh)/(atoms%rmsh(:, n)*atoms%rmsh(:, n))
133 :
134 72246 : DO js = 1, jspV
135 : chlh(1:atoms%jri(n), lh, js) = den_mt(1:atoms%jri(n), lh, js) / &
136 37164800 : (atoms%rmsh(1:atoms%jri(n), n)*atoms%rmsh(1:atoms%jri(n), n))
137 :
138 : !Necessary gradients
139 71554 : IF (dograds) THEN
140 : !Colinear case only needs radial derivatives of chlh
141 15820 : CALL grdchlh(atoms%dx(n), chlh(1:atoms%jri(n), lh, js), chlhdr(1:, lh, js), chlhdrr(1:, lh,js),atoms%rmsh(:,n))
142 41963 : IF (any(noco%l_unrestrictMT)) THEN
143 : !Noco case also needs radial derivatives of mm
144 0 : CALL grdchlh(atoms%dx(n), mm(:atoms%jri(n),lh), drm(:,lh), drrm(:,lh), atoms%rmsh(:, n))
145 : END IF
146 : END IF
147 : END DO ! js
148 : END DO ! lh
149 :
150 : !The following Loop maps chlh on the k-Grid using the lattice harmonics ylh
151 : !$OMP parallel do default( NONE ) &
152 : !$OMP SHARED(atoms,sphhar,noco,n,nsp,jspV,nd,ylh,chlh,dograds,chlhdr,chlhdrr,mm,drm,drrm)&
153 : !$OMP SHARED(ylht,ylhf,ylhtt,ylhff,ylhtf,jspins,thet,grad,ch,ch_calc) &
154 : !$OMP private(kt,ch_tmp,js,lh,k,chdr,chdt,chdf,chdrr,chdtt,chdff,chdtf,chdrt,chdrf) &
155 692 : !$OMP private(chlhtot,chlhdrtot,chlhdrrtot)
156 : DO jr = 1, atoms%jri(n)
157 : kt = (jr-1)*nsp
158 : ! charge density (on extended grid for all jr)
159 : ! following are at points on jr-th sphere.
160 : ch_tmp(:, :) = 0.0
161 : ! generate the densities on an angular mesh (ch_tmp is needed in mkgylm call later on)
162 : DO js = 1, jspV
163 : DO lh = 0, sphhar%nlh(nd)
164 : DO k = 1, nsp
165 : ch_tmp(k, js) = ch_tmp(k, js) + ylh(k, lh, nd)*chlh(jr, lh, js)
166 : ENDDO
167 : ENDDO
168 : ENDDO
169 :
170 : !Initialize derivatives of ch on grid if needed.
171 : IF (dograds) THEN
172 : chdr(:, :) = 0.0 ! d(ch)/dr
173 : chdt(:, :) = 0.0 ! d(ch)/dtheta
174 : chdf(:, :) = 0.0 ! d(ch)/dfai
175 : chdrr(:, :) = 0.0 ! dd(ch)/drr
176 : chdtt(:, :) = 0.0 ! dd(ch)/dtt
177 : chdff(:, :) = 0.0 ! dd(ch)/dff
178 : chdtf(:, :) = 0.0 ! dd(ch)/dtf
179 : chdrt(:, :) = 0.0 ! d(d(ch)/dr)dt
180 : chdrf(:, :) = 0.0 ! d(d(ch)/dr)df
181 : ! generate the derivatives on an angular mesh
182 : DO js = 1, jspV
183 : DO lh = 0, sphhar%nlh(nd)
184 :
185 : !The following snippet maps chlh and its radial derivatives on a colinear system
186 : !using mm and its radial derivatives.
187 : IF (any(noco%l_unrestrictMT)) THEN
188 : IF (js.EQ.1) THEN
189 : chlhtot(lh)=0.5*(chlh(jr, lh, 1)+chlh(jr, lh, 2))
190 : chlhdrtot(lh)=0.5*(chlhdr(jr, lh, 1)+chlhdr(jr, lh, 2))
191 : chlhdrrtot(lh)=0.5*(chlhdrr(jr, lh, 1)+chlhdrr(jr, lh, 2))
192 : chlh(jr,lh,js)=chlhtot(lh)+mm(jr,lh)
193 : chlhdr(jr, lh, js)=chlhdrtot(lh)+drm(jr,lh)
194 : chlhdrr(jr, lh, js)=chlhdrrtot(lh)+drrm(jr,lh)
195 : ELSE IF (js.EQ.2) THEN
196 : chlh(jr,lh,js)=chlhtot(lh)-mm(jr,lh)
197 : chlhdr(jr, lh, js)=chlhdrtot(lh)-drm(jr,lh)
198 : chlhdrr(jr, lh, js)=chlhdrrtot(lh)-drrm(jr,lh)
199 : ELSE
200 : chlh(jr,lh,js)=0
201 : chlhdr(jr, lh, js)=0
202 : chlhdrr(jr, lh, js)=0
203 : END IF
204 : END IF
205 :
206 : !The following loop brings chlhdr and chlhdrr on the k-grid.
207 : DO k = 1, nsp
208 : chdr(k, js) = chdr(k, js) + ylh(k, lh, nd)*chlhdr(jr, lh, js)
209 : chdrr(k, js) = chdrr(k, js) + ylh(k, lh, nd)*chlhdrr(jr, lh, js)
210 : ENDDO
211 :
212 : !This loop calculates the other derviatives of ch (Angular terms) on the k-grid
213 : !by using the lattice harmonics derivatives and chlh with its derivatives.
214 : DO k = 1, nsp
215 : chdrt(k, js) = chdrt(k, js) + ylht(k, lh, nd)*chlhdr(jr, lh, js)
216 : chdrf(k, js) = chdrf(k, js) + ylhf(k, lh, nd)*chlhdr(jr, lh, js)
217 : chdt(k, js) = chdt(k, js) + ylht(k, lh, nd)*chlh(jr, lh, js)
218 : chdf(k, js) = chdf(k, js) + ylhf(k, lh, nd)*chlh(jr, lh, js)
219 : chdtt(k, js) = chdtt(k, js) + ylhtt(k, lh, nd)*chlh(jr, lh, js)
220 : chdff(k, js) = chdff(k, js) + ylhff(k, lh, nd)*chlh(jr, lh, js)
221 : chdtf(k, js) = chdtf(k, js) + ylhtf(k, lh, nd)*chlh(jr, lh, js)
222 : ENDDO
223 : ENDDO ! lh
224 : ENDDO ! js
225 : !Rotation to local if needed (Indicated by rotch)
226 : !Makegradients
227 : !IF(jspins>2) CALL mkgylm(2, atoms%rmsh(jr, n), thet, nsp, &
228 : ! ch_tmp, chdr, chdt, chdf, chdrr, chdtt, chdff, chdtf, chdrt, chdrf, grad, kt)
229 : !IF(jspins.LE.2)
230 : CALL mkgylm(jspins, atoms%rmsh(jr, n), thet, nsp, &
231 : ch_tmp, chdr, chdt, chdf, chdrr, chdtt, chdff, chdtf, chdrt, chdrf, grad, kt)
232 : END IF
233 : !Set charge to minimum value
234 : IF (PRESENT(ch)) THEN
235 : WHERE (ABS(ch_tmp(:nsp,:)) < d_15) ch_tmp(:nsp,:) = d_15
236 : ch_calc(kt + 1:kt + nsp, :) = ch_tmp(:nsp, :)
237 : ENDIF
238 : END DO
239 : !$OMP END PARALLEL DO
240 :
241 :
242 692 : IF (PRESENT(ch)) THEN
243 : !Rotation to local if needed (Indicated by rotch)
244 1947 : IF (rotch.AND.any(noco%l_unrestrictMT)) THEN
245 2247301 : DO jr = 1,nsp*atoms%jri(n)
246 2247284 : rho_11 = ch_calc(jr,1)
247 2247284 : rho_22 = ch_calc(jr,2)
248 2247284 : rho_21r = ch_calc(jr,3)
249 2247284 : rho_21i = ch_calc(jr,4)
250 2247284 : mx = 2*rho_21r
251 2247284 : my = -2*rho_21i
252 2247284 : mz = (rho_11-rho_22)
253 2247284 : magmom = SQRT(mx**2 + my**2 + mz**2)
254 2247284 : rhotot = rho_11 + rho_22
255 2247284 : rho_up = (rhotot + magmom)/2
256 2247284 : rho_down= (rhotot - magmom)/2
257 2247284 : ch(jr,1) = rho_up
258 2247301 : ch(jr,2) = rho_down
259 : END DO
260 : ELSE
261 149000332 : ch(:nsp*atoms%jri(n),1:jspins)=ch_calc(:nsp*atoms%jri(n),1:jspins)
262 :
263 : END IF
264 : END IF
265 692 : call timestop("mt_to_grid")
266 692 : END SUBROUTINE mt_to_grid
267 :
268 2522 : SUBROUTINE mt_from_grid(atoms, sym, sphhar, n, jspins, v_in, vr)
269 : IMPLICIT NONE
270 : TYPE(t_atoms), INTENT(IN) :: atoms
271 : TYPE(t_sym), INTENT(IN) :: sym
272 : TYPE(t_sphhar), INTENT(IN):: sphhar
273 : INTEGER, INTENT(IN) :: jspins, n
274 : REAL, INTENT(IN) :: v_in(:, :)
275 : REAL, INTENT(INOUT) :: vr(:, 0:, :)
276 :
277 2522 : REAL :: vpot(atoms%nsp()), vlh
278 : INTEGER :: js, kt, lh, jr, nd, nsp
279 :
280 2522 : call timestart("mt_from_grid")
281 :
282 2522 : nsp = atoms%nsp()
283 2522 : nd = sym%ntypsy(atoms%firstAtom(n))
284 :
285 6215 : DO js = 1, jspins
286 : !$OMP parallel do default( NONE ) &
287 : !$OMP SHARED(atoms,sphhar,n,v_in,nsp,wt,nd,ylh,vr,js)&
288 6215 : !$OMP private(vpot,kt,lh,vlh)
289 : DO jr = 1, atoms%jri(n)
290 : kt = (jr-1)*nsp
291 : vpot = v_in(kt + 1:kt + nsp, js)*wt(:)! multiplicate v_in with the weights of the k-points
292 :
293 : DO lh = 0, sphhar%nlh(nd)
294 : !
295 : ! ---> determine the corresponding potential number
296 : !c through gauss integration
297 : !
298 : vlh = dot_PRODUCT(vpot(:), ylh(:nsp, lh, nd))
299 : vr(jr, lh, js) = vr(jr, lh, js) + vlh
300 : ENDDO ! lh
301 : ENDDO ! jr
302 : !$OMP END PARALLEL DO
303 : ENDDO
304 2522 : call timestop("mt_from_grid")
305 2522 : END SUBROUTINE mt_from_grid
306 :
307 732 : SUBROUTINE finish_mt_grid()
308 : implicit NONE
309 732 : call timestart("finish_mt_grid")
310 732 : DEALLOCATE (ylh, wt, rx, thet, phi)
311 732 : IF (ALLOCATED(ylht)) DEALLOCATE (ylht, ylhtt, ylhf, ylhff, ylhtf)
312 732 : call timestop("finish_mt_grid")
313 732 : END SUBROUTINE finish_mt_grid
314 :
315 : END MODULE m_mt_tofrom_grid
|
