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 :
7 : MODULE m_rotate_int_den_tofrom_local
8 : USE m_juDFT
9 : USE m_fft2d
10 : USE m_fft3d
11 : USE m_types
12 :
13 : IMPLICIT NONE
14 :
15 : CONTAINS
16 :
17 192 : SUBROUTINE rotate_int_den_to_local(sym, stars, atoms, sphhar, vacuum, cell, &
18 : input, noco, den)
19 :
20 : !--------------------------------------------------------------------------
21 : ! This subroutine calculates the spin-up and -down density in the intersti-
22 : ! tial region i.e. it takes the non-collinear density and rotates it into
23 : ! a local spin-frame, making it spin-diagonal.
24 : !
25 : ! The rotated density is needed to calculate the potential-energy integrals
26 : ! in vgen_xcpot. For accuracy reasons, the magnetisation for the potential
27 : ! itself is regenerated from the unrotated densities.
28 : !
29 : ! In addition this routine stores the angles used in the rotation. These
30 : ! angles are later needed to rotate the up- and down-potentials back to the
31 : ! global frame. DW 2018
32 : !
33 : ! Based on rhodirgen by
34 : ! Philipp Kurz 99/11/01
35 : !--------------------------------------------------------------------------
36 :
37 : !-------Important variables:-----------------------------------------------
38 : ! ifft3: size of the 3d real space mesh
39 : ! ifft2: size of the 2d real space mesh
40 : ! ris: first components of the density matrix
41 : ! later interstitial spin-up and -down density and direction of mag-
42 : ! netic field (theta and phi) all stored on real space mesh
43 : !--------------------------------------------------------------------------
44 :
45 : USE m_constants
46 : USE m_polangle
47 :
48 : TYPE(t_noco), INTENT(IN) :: noco
49 :
50 : TYPE(t_input), INTENT(IN) :: input
51 : TYPE(t_vacuum), INTENT(IN) :: vacuum
52 : TYPE(t_sym), INTENT(IN) :: sym
53 : TYPE(t_stars), INTENT(IN) :: stars
54 : TYPE(t_cell), INTENT(IN) :: cell
55 : TYPE(t_sphhar), INTENT(IN) :: sphhar
56 : TYPE(t_atoms), INTENT(IN) :: atoms
57 : TYPE(t_potden), INTENT(INOUT) :: den
58 :
59 : INTEGER :: iden, jspin, ivac, ifft2, ifft3
60 : INTEGER :: imz, ityp, iri, ilh, imesh, iq2, iq3
61 : REAL :: rho_11, rho_22, rho_21r, rho_21i
62 : REAL :: mx, my, mz, magmom, vz_r, vz_i, rziw
63 : REAL :: rhotot, rho_up, rho_down, theta, phi
64 : REAL :: eps=1E-20
65 :
66 192 : REAL, ALLOCATABLE :: rz(:,:,:), rvacxy(:,:,:,:), ris(:,:)
67 192 : REAL, ALLOCATABLE :: fftwork(:)
68 :
69 : ! Initialize arrays for the density matrix:
70 :
71 192 : ifft3 = 27*stars%mx1*stars%mx2*stars%mx3
72 192 : IF (input%film) THEN
73 0 : ifft2 = 9*stars%mx1*stars%mx2
74 : ELSE
75 : ifft2=0
76 : END IF
77 :
78 192 : IF (ALLOCATED(den%phi_pw)) THEN
79 0 : DEALLOCATE(den%phi_pw)!,den%phi_vacz,den%phi_vacxy)
80 0 : DEALLOCATE(den%theta_pw)!,den%theta_vacz,den%theta_vacxy)
81 0 : DEALLOCATE(den%theta_vac,den%phi_vac)
82 : END IF
83 :
84 768 : ALLOCATE(den%phi_pw(ifft3),den%theta_pw(ifft3))
85 : !ALLOCATE(den%phi_vacz(vacuum%nmzd,2),den%theta_vacz(vacuum%nmzd,2))
86 1344 : ALLOCATE(den%phi_vac(ifft2,vacuum%nmzd,2),den%theta_vac(ifft2,vacuum%nmzd,2))
87 : !ALLOCATE(den%phi_vacxy(ifft2,vacuum%nmzxyd,2),den%theta_vacxy(ifft2,vacuum%nmzxyd,2))
88 :
89 768 : ALLOCATE (ris(ifft3,4),fftwork(ifft3))
90 :
91 : ! Interstitial part:
92 :
93 : ! Fourier transform the diagonal part of the density matrix (den%pw)
94 : ! to real space (ris):
95 576 : DO iden = 1,2
96 576 : CALL fft3d(ris(:,iden),fftwork,den%pw(:,iden),stars,+1)
97 : END DO
98 :
99 : ! Fourier transform the off-diagonal part of the density matrix:
100 192 : CALL fft3d(ris(:,3),ris(:,4),den%pw(:,3),stars,+1)
101 :
102 : ! Calculate the charge and magnetization densities on the real space mesh:
103 7185054 : DO imesh = 1,ifft3
104 7184862 : rho_11 = ris(imesh,1)
105 7184862 : rho_22 = ris(imesh,2)
106 7184862 : rho_21r = ris(imesh,3)
107 7184862 : rho_21i = ris(imesh,4)
108 7184862 : mx = 2*rho_21r
109 7184862 : my = -2*rho_21i ! TODO: This is a magic minus.
110 7184862 : mz = rho_11 - rho_22
111 7184862 : magmom = SQRT(mx**2 + my**2 + mz**2)
112 7184862 : rhotot = rho_11 + rho_22
113 7184862 : rho_up = (rhotot + magmom)/2
114 7184862 : rho_down = (rhotot - magmom)/2
115 :
116 7184862 : CALL pol_angle(mx,my,mz,theta,phi)
117 :
118 7184862 : ris(imesh,1) = rho_up
119 7184862 : ris(imesh,2) = rho_down
120 7184862 : den%theta_pw(imesh) = theta
121 7185054 : den%phi_pw(imesh) = phi
122 : END DO
123 :
124 : ! Fourier transform the matrix potential back to reciprocal space:
125 576 : DO jspin = 1, input%jspins
126 14370108 : fftwork=0.0
127 576 : CALL fft3d(ris(:,jspin),fftwork,den%pw(:,jspin),stars,-1)
128 : END DO
129 :
130 192 : IF (.NOT.input%film) RETURN
131 :
132 : ! Now the vacuum part starts:
133 :
134 0 : ALLOCATE(rvacxy(ifft2,vacuum%nmzxyd,2,4))
135 0 : ALLOCATE (rz(vacuum%nmzd,2,2))
136 :
137 : ! Fourier transform the diagonal part of the density matrix (den%vacz and
138 : ! den%vacxy) to real space (rvacxy):
139 0 : DO iden = 1,2
140 0 : DO ivac = 1,vacuum%nvac
141 0 : DO imz = 1,vacuum%nmzxyd
142 0 : rziw = 0.0
143 : CALL fft2d(stars,rvacxy(:,imz,ivac,iden),fftwork,&
144 : den%vac(imz,:,ivac,iden),&
145 0 : 1)
146 : END DO
147 : END DO
148 : END DO
149 :
150 : ! Fourier transform the off-diagonal part of the density matrix:
151 0 : DO ivac = 1,vacuum%nvac
152 0 : DO imz = 1,vacuum%nmzxyd
153 0 : rziw = 0.0
154 0 : vz_r = REAL(den%vac(imz,1,ivac,3))
155 0 : vz_i = AIMAG(den%vac(imz,1,ivac,3))
156 : CALL fft2d(stars,rvacxy(:,imz,ivac,3),rvacxy(:,imz,ivac,4),&
157 0 : den%vac(imz,:,ivac,3),1)
158 :
159 : END DO
160 : END DO
161 :
162 : ! Calculate the charge and magnetization densities on the real space mesh:
163 0 : DO ivac = 1,vacuum%nvac
164 0 : DO imz = 1,vacuum%nmzxyd
165 0 : DO imesh = 1,ifft2
166 0 : rho_11 = rvacxy(imesh,imz,ivac,1)
167 0 : rho_22 = rvacxy(imesh,imz,ivac,2)
168 0 : rho_21r = rvacxy(imesh,imz,ivac,3)
169 0 : rho_21i = rvacxy(imesh,imz,ivac,4)
170 0 : mx = 2*rho_21r
171 0 : my = -2*rho_21i
172 0 : mz = rho_11 - rho_22
173 0 : magmom = SQRT(mx**2 + my**2 + mz**2)
174 0 : rhotot = rho_11 + rho_22
175 0 : rho_up = (rhotot + magmom)/2
176 0 : rho_down = (rhotot - magmom)/2
177 :
178 0 : CALL pol_angle(mx,my,mz,theta,phi)
179 :
180 0 : rvacxy(imesh,imz,ivac,1) = rho_up
181 0 : rvacxy(imesh,imz,ivac,2) = rho_down
182 0 : den%theta_vac(imesh,imz,ivac) = theta
183 0 : den%phi_vac(imesh,imz,ivac) = phi
184 : END DO
185 : END DO
186 :
187 0 : DO imz = vacuum%nmzxyd+1,vacuum%nmzd
188 0 : rho_11 = REAL(den%vac(imz,1,ivac,1))
189 0 : rho_22 = REAL(den%vac(imz,1,ivac,2))
190 0 : rho_21r = REAL(den%vac(imz,1,ivac,3))
191 0 : rho_21i = AIMAG(den%vac(imz,1,ivac,3))
192 0 : mx = 2*rho_21r
193 0 : my = -2*rho_21i
194 0 : mz = rho_11 - rho_22
195 0 : magmom = SQRT(mx**2 + my**2 + mz**2)
196 0 : rhotot = rho_11 + rho_22
197 0 : rho_up = (rhotot + magmom)/2
198 0 : rho_down = (rhotot - magmom)/2
199 :
200 0 : CALL pol_angle(mx,my,mz,theta,phi)
201 :
202 0 : den%vac(imz,1,ivac,1) = rho_up
203 0 : den%vac(imz,1,ivac,2) = rho_down
204 0 : den%theta_vac(1,imz,ivac) = theta
205 0 : den%phi_vac(1,imz,ivac) = phi
206 : END DO
207 : END DO
208 :
209 : ! Fourier transform the matrix potential back to reciprocal space:
210 0 : DO jspin = 1,input%jspins
211 0 : DO ivac = 1,vacuum%nvac
212 0 : DO imz = 1,vacuum%nmzxyd
213 0 : fftwork=0.0
214 : CALL fft2d(stars,rvacxy(:,imz,ivac,jspin),fftwork,&
215 : den%vac(imz,:,ivac,jspin),&
216 0 : -1)
217 :
218 : END DO
219 : END DO
220 : END DO
221 :
222 : RETURN
223 :
224 192 : END SUBROUTINE rotate_int_den_to_local
225 :
226 192 : SUBROUTINE rotate_int_den_from_local(stars,atoms,vacuum,sym,input,den,vTot)
227 :
228 : !--------------------------------------------------------------------------
229 : ! This subroutine prepares the spin dependent 2x2 matrix potential for the
230 : ! Hamiltonian setup. This is done in 4 steps.
231 : !
232 : ! i) The spin up and down potential and the angles of the magentic field,
233 : ! theta and phi, are reloaded from den.
234 : ! ii) The spin up and down potential is Fourier transformed to real space
235 : ! (theta and phi are also stored on the real space grid).
236 : ! iii) The four components of the matrix potential are calculated on the
237 : ! real space mesh.
238 : ! iv) The matrix potential is Fourier transformed, stored in terms of
239 : ! stars and written to vTot%pw(_w).
240 : !
241 : ! Philipp Kurz 99/11/01
242 : !--------------------------------------------------------------------------
243 :
244 : !-------Important variables:-----------------------------------------------
245 : ! ifft3: size of the 3d real space mesh
246 : ! ifft2: size of the 2d real space mesh
247 : ! vis: first interstitial spin up and down potential and angles of magnetic
248 : ! field (theta and phi)
249 : ! later four components of matrix potential all stored in real space
250 : !--------------------------------------------------------------------------
251 :
252 : !
253 : TYPE(t_input), INTENT(IN) :: input
254 : TYPE(t_vacuum), INTENT(IN) :: vacuum
255 : TYPE(t_sym), INTENT(IN) :: sym
256 : TYPE(t_stars), INTENT(IN) :: stars
257 : TYPE(t_atoms), INTENT(IN) :: atoms
258 : TYPE(t_potden), INTENT(IN) :: den
259 : TYPE(t_potden), INTENT(INOUT) :: vTot
260 :
261 : INTEGER :: imeshpt, ipot, jspin, ig2, ig3, ivac
262 : INTEGER :: ifft2, ifft3, imz, iter, i
263 : REAL :: vup, vdown, veff, beff, vziw, theta, phi
264 :
265 192 : REAL, ALLOCATABLE :: vvacxy(:,:,:,:), vis(:,:), vis2(:,:)
266 : REAL, ALLOCATABLE :: fftwork(:)
267 :
268 : ! Initialize arrays for the potential matrix:
269 :
270 192 : ifft3 = 27*stars%mx1*stars%mx2*stars%mx3
271 192 : IF (ifft3.NE.SIZE(den%theta_pw)) CALL judft_error("Wrong size of angles")
272 192 : ifft2 = SIZE(den%phi_vac,1)
273 :
274 1152 : ALLOCATE ( vis(ifft3,4),fftwork(ifft3),vis2(ifft3,4))
275 :
276 : ! Interstitial part:
277 :
278 : ! Fourier transform the diagonal part of the potential matrix (vTot%pw)
279 : ! to real space (vis):
280 576 : DO jspin = 1,input%jspins
281 576 : CALL fft3d(vis(:,jspin),fftwork, vTot%pw(:,jspin), stars,+1)
282 : END DO
283 :
284 : ! Calculate the four components of the matrix potential on the real space
285 : ! mesh:
286 7185054 : DO imeshpt = 1, ifft3
287 7184862 : vup = vis(imeshpt,1)
288 7184862 : vdown = vis(imeshpt,2)
289 7184862 : theta = den%theta_pw(imeshpt)
290 7184862 : phi = den%phi_pw(imeshpt)
291 :
292 7184862 : veff = (vup + vdown)/2.0
293 7184862 : beff = (vup - vdown)/2.0
294 :
295 7184862 : vis(imeshpt,1) = veff + beff*COS(theta) ! V_(1,1) [V+B_z]
296 7184862 : vis(imeshpt,2) = veff - beff*COS(theta) ! V_(2,2) [V-B_z]
297 7184862 : vis(imeshpt,3) = beff*SIN(theta)*COS(phi) ! Re(V_(2,1)) [B_x]
298 7184862 : vis(imeshpt,4) = beff*SIN(theta)*SIN(phi) ! Im(V_(2,1)) [B_y]
299 :
300 35924502 : DO ipot = 1,4
301 35924310 : vis2(imeshpt,ipot) = vis(imeshpt,ipot) * stars%ufft(imeshpt-1)
302 : END DO
303 : END DO
304 :
305 : ! Fourier transform the matrix potential back to reciprocal space:
306 576 : DO ipot = 1,2
307 14370108 : fftwork=0.0
308 384 : CALL fft3d(vis(:,ipot),fftwork, vTot%pw(1,ipot), stars,-1)
309 14370108 : fftwork=0.0
310 576 : CALL fft3d(vis2(:,ipot),fftwork, vTot%pw_w(1,ipot), stars,-1)
311 : END DO
312 :
313 192 : CALL fft3d(vis(:,3),vis(:,4), vTot%pw(1,3), stars,-1)
314 192 : CALL fft3d(vis2(:,3),vis2(:,4), vTot%pw_w(1,3), stars,-1)
315 :
316 192 : IF (.NOT. input%film) RETURN
317 :
318 : ! Now the vacuum part starts:
319 :
320 0 : ALLOCATE(vvacxy(ifft2,vacuum%nmzxyd,2,4))
321 :
322 : ! Fourier transform the up and down potentials (vTot%vacz and vTot%vacxy)
323 : ! to real space (vvacxy):
324 0 : DO jspin = 1,input%jspins
325 0 : DO ivac = 1,vacuum%nvac
326 0 : DO imz = 1,vacuum%nmzxyd
327 0 : vziw = 0.0
328 : !
329 : CALL fft2d(stars, vvacxy(:,imz,ivac,jspin),fftwork,&
330 0 : vTot%vac(imz,:,ivac,jspin), 1)
331 :
332 : END DO
333 : END DO
334 : END DO
335 :
336 : ! Calculate the four components of the matrix potential in real space:
337 0 : DO ivac = 1,vacuum%nvac
338 0 : DO imz = 1,vacuum%nmzxyd
339 0 : DO imeshpt = 1,ifft2
340 0 : vup = vvacxy(imeshpt,imz,ivac,1)
341 0 : vdown = vvacxy(imeshpt,imz,ivac,2)
342 0 : theta = den%theta_vac(imeshpt,imz,ivac)
343 0 : phi = den%phi_vac(imeshpt,imz,ivac)
344 :
345 0 : veff = (vup + vdown)/2.0
346 0 : beff = (vup - vdown)/2.0
347 0 : vvacxy(imeshpt,imz,ivac,1) = veff + beff*COS(theta)
348 0 : vvacxy(imeshpt,imz,ivac,2) = veff - beff*COS(theta)
349 0 : vvacxy(imeshpt,imz,ivac,3) = beff*SIN(theta)*COS(phi)
350 0 : vvacxy(imeshpt,imz,ivac,4) = beff*SIN(theta)*SIN(phi)
351 : END DO
352 : END DO
353 :
354 0 : DO imz = vacuum%nmzxyd+1,vacuum%nmzd
355 0 : vup = REAL(vTot%vac(imz,1,ivac,1))
356 0 : vdown = REAL(vTot%vac(imz,1,ivac,2))
357 0 : theta = den%theta_vac(1,imz,ivac)
358 0 : phi = den%phi_vac(1,imz,ivac)
359 0 : veff = (vup + vdown)/2.0
360 0 : beff = (vup - vdown)/2.0
361 0 : vTot%vac(imz,1,ivac,1) = veff + beff*COS(theta)
362 0 : vTot%vac(imz,1,ivac,2) = veff - beff*COS(theta)
363 0 : vTot%vac(imz,1,ivac,3) = beff*SIN(theta)*COS(phi)+ImagUnit*beff*SIN(theta)*SIN(phi)
364 : END DO
365 : END DO
366 :
367 : ! Fourier transform the matrix potential back to reciprocal space:
368 0 : DO ipot = 1,2
369 0 : DO ivac = 1,vacuum%nvac
370 0 : DO imz = 1,vacuum%nmzxyd
371 0 : fftwork=0.0
372 : !
373 : CALL fft2d(stars, vvacxy(:,imz,ivac,ipot),fftwork,&
374 0 : vTot%vac(imz,:,ivac,ipot),-1)
375 :
376 : END DO
377 : END DO
378 : END DO
379 :
380 0 : DO ivac = 1,vacuum%nvac
381 0 : DO imz = 1,vacuum%nmzxyd
382 0 : fftwork=0.0
383 : !
384 : CALL fft2d(stars, vvacxy(:,imz,ivac,3),vvacxy(:,imz,ivac,4),&
385 0 : vTot%vac(imz,:,ivac,3),-1)
386 :
387 : END DO
388 : END DO
389 :
390 : RETURN
391 :
392 192 : END SUBROUTINE rotate_int_den_from_local
393 :
394 : END MODULE m_rotate_int_den_tofrom_local
|
