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_pw_tofrom_grid
7 : USE m_types
8 : USE m_types_fftGrid
9 : PRIVATE
10 : REAL,PARAMETER:: d_15=1.e-15
11 :
12 : TYPE(t_fftgrid) :: fftgrid
13 : INTEGER :: griddim
14 : REAL :: gmax
15 : REAL :: gmaxz
16 :
17 : PUBLIC :: init_pw_grid, pw_to_grid, pw_from_grid, finish_pw_grid
18 : CONTAINS
19 339 : SUBROUTINE init_pw_grid(stars,sym,cell,xcpot)
20 : USE m_types
21 : use m_types_xcpot
22 : IMPLICIT NONE
23 : TYPE(t_stars),INTENT(IN) :: stars
24 : TYPE(t_sym),INTENT(IN) :: sym
25 : TYPE(t_cell),INTENT(IN) :: cell
26 : CLASS(t_xcpot),INTENT(IN),OPTIONAL :: xcpot
27 :
28 : !---> set up pointer for backtransformation of from g-vector in
29 : ! positive domain of xc density fftbox into stars.
30 : ! also the x,y,z components of the g-vectors are set up to calculate
31 : ! derivatives.
32 : ! in principle this can also be done in main program once.
33 : ! it is done here to save memory.
34 :
35 339 : if (present(xcpot)) THEN
36 337 : gmax=xcpot%gmaxxc
37 337 : gmaxz=stars%gmaxz
38 337 : if (xcpot%needs_grad().and.gmax>0.0) then
39 276 : if (gmaxz>0.0) then
40 0 : call fftgrid%init(cell,sym,gmax,gmaxz)
41 : else
42 276 : call fftgrid%init(cell,sym,gmax)
43 : end if
44 : else
45 244 : call fftgrid%init((/3*stars%mx1,3*stars%mx2,3*stars%mx3/))
46 61 : gmax=stars%gmax
47 : endif
48 : else
49 2 : gmax=stars%gmax
50 2 : gmaxz=stars%gmaxz
51 2 : if (gmaxz>0.0) then
52 0 : call fftgrid%init(cell,sym,gmax,gmaxz)
53 : else
54 2 : call fftgrid%init(cell,sym,gmax)
55 : end if
56 : endif
57 339 : griddim=size(fftgrid%grid)
58 :
59 339 : END SUBROUTINE init_pw_grid
60 :
61 702 : SUBROUTINE pw_to_grid(dograds,jspins,l_noco,stars,cell,den_pw,grad,xcpot,rho,rhoim)
62 : !.....------------------------------------------------------------------
63 : !-------> abbreviations
64 : !
65 : ! ph_wrk: work array containing phase * g_x,gy......
66 : ! den%pw: charge density stored as stars
67 : ! rho : charge density stored in real space
68 : ! v_xc : exchange-correlation potential in real space
69 : ! exc : exchange-correlation energy density in real space
70 : ! kxc1d : dimension of the charge density fft box in the pos. domain
71 : ! kxc2d : defined in dimens.f program (subroutine apws).1,2,3 indic
72 : ! kxc3d ; a_1, a_2, a_3 directions.
73 : ! kq(i) : i=1,2,3 actual length of the fft-box for which fft is done
74 : ! nstr : number of members (arms) of reciprocal lattice (g) vector
75 : ! of each star.
76 : ! nxc3_fft: number of stars in the charge density fft-box
77 : ! ng3 : number of 3 dim. stars in the charge density sphere define
78 : ! by gmax
79 : ! kmxxc_fft: number of g-vectors forming the nxc3_fft stars in the
80 : ! charge density or xc-density sphere
81 : ! kimax : number of g-vectors forming the ng3 stars in the gmax-sphe
82 : ! igfft : pointer from the g-sphere (stored as stars) to fft-grid
83 : ! and from fft-grid to g-sphere (stored as stars)
84 : ! pgfft : contains the phases of the g-vectors of sph.
85 : ! isn : isn = +1, fft transform for g-space to r-space
86 : ! isn = -1, vice versa
87 : ! l_rdm : if true, calculate noco gradients in local frame of density matrix
88 : ! improves convergence in case of GGA + noco (best with large G_max_xc)
89 : !
90 : !-------------------------------------------------------------------
91 : USE m_grdrsis
92 : USE m_polangle
93 : USE m_mkgxyz3
94 : USE m_types
95 : USE m_constants
96 : IMPLICIT NONE
97 :
98 : LOGICAL,INTENT(IN) :: dograds
99 : INTEGER,INTENT(IN) :: jspins
100 : LOGICAL,INTENT(IN) :: l_noco
101 : TYPE(t_stars),INTENT(IN) :: stars
102 : TYPE(t_cell),INTENT(IN) :: cell
103 : COMPLEX,INTENT(IN) :: den_pw(:,:)
104 : TYPE(t_gradients),INTENT(OUT) :: grad
105 : CLASS(t_xcpot), INTENT(IN),OPTIONAL :: xcpot
106 : REAL,ALLOCATABLE,INTENT(OUT),OPTIONAL :: rho(:,:),rhoim(:,:)
107 :
108 :
109 : INTEGER :: js,i,idm,ig,ndm,jdm,j
110 : REAL :: rhotot,mmx,mmy,mmz,theta,phi,fd(3),sd(3), rhoCutGrad
111 : COMPLEX :: rho21
112 : ! .. Local Arrays ..
113 351 : COMPLEX, ALLOCATABLE :: cqpw(:,:),ph_wrk(:)
114 351 : REAL, ALLOCATABLE :: bf3(:)
115 351 : REAL, ALLOCATABLE :: rhd1(:,:,:),rhd2(:,:,:)
116 351 : REAL, ALLOCATABLE :: mx(:),my(:)
117 351 : REAL, ALLOCATABLE :: magmom(:),dmagmom(:,:),ddmagmom(:,:,:)
118 351 : REAL, ALLOCATABLE :: rhodiag(:,:),der(:,:,:),dder(:,:,:,:)
119 351 : REAL, ALLOCATABLE :: sinsqu(:),cossqu(:),sincos(:),rhdd(:,:,:,:)
120 351 : COMPLEX, ALLOCATABLE :: exi(:)
121 :
122 : LOGICAL, PARAMETER :: l_rdm=.true.
123 :
124 : ! Allocate arrays
125 351 : ALLOCATE( bf3(0:griddim-1))
126 351 : IF (dograds) THEN
127 560 : IF (PRESENT(rho)) ALLOCATE(rho(0:griddim-1,jspins))
128 284 : IF (PRESENT(rhoim)) ALLOCATE(rhoim(0:griddim-1,jspins))
129 852 : ALLOCATE( ph_wrk(0:griddim-1),rhd1(0:griddim-1,jspins,3))
130 852 : ALLOCATE( rhd2(0:griddim-1,jspins,6) )
131 : ELSE
132 134 : IF (PRESENT(rho)) ALLOCATE(rho(0:griddim-1,jspins))
133 67 : IF (PRESENT(rhoim)) ALLOCATE(rhoim(0:griddim-1,jspins))
134 : ENDIF
135 351 : IF (l_noco) THEN
136 96 : IF (dograds) THEN
137 67 : ALLOCATE( mx(0:griddim-1),my(0:griddim-1),magmom(0:griddim-1))
138 : IF (l_rdm) THEN
139 134 : ALLOCATE( rhodiag(0:griddim-1,jspins),der(0:griddim-1,3,4),dder(0:griddim-1,3,3,4),rhdd(0:griddim-1,2,3,3) )
140 67 : ALLOCATE( sinsqu(0:griddim-1),cossqu(0:griddim-1),sincos(0:griddim-1),exi(0:griddim-1) )
141 : ELSE
142 : ALLOCATE( dmagmom(0:griddim-1,3),ddmagmom(0:griddim-1,3,3) )
143 : ENDIF
144 : ELSE
145 29 : ALLOCATE( mx(0:griddim-1),my(0:griddim-1),magmom(0:griddim-1))
146 : ENDIF
147 : END IF
148 : ! if (gmaxz>0.0) write(*,*) "Present in pw_to_grid"
149 351 : IF (PRESENT(rho)) THEN
150 : !Put den_pw on grid and store into rho(:,1:2)
151 885 : DO js=1,jspins
152 542 : if (gmaxz>0.0) then
153 0 : call fftgrid%putFieldOnGrid(stars,den_pw(:,js),cell,gmax,gmaxz)
154 : else
155 542 : call fftgrid%putFieldOnGrid(stars,den_pw(:,js),cell,gmax)
156 : end if
157 542 : call fftgrid%perform_fft(forward=.false.)
158 : ! TODO: grid is technically still complex right? The REAL cast happens here:
159 : !rho(0:,js)=fftgrid%grid
160 10085863 : rho(0:,js) = REAL(fftgrid%grid)
161 885 : IF (PRESENT(rhoim)) rhoim(0:,js) = AIMAG(fftgrid%grid)
162 : END DO
163 :
164 343 : IF (l_noco) THEN
165 : ! Get mx,my on real space grid and recalculate rho and magmom
166 96 : if (gmaxz>0.0) then
167 0 : call fftgrid%putFieldOnGrid(stars,den_pw(:,3),cell,gmax,gmaxz)
168 : else
169 96 : call fftgrid%putFieldOnGrid(stars,den_pw(:,3),cell,gmax)
170 : end if
171 96 : call fftgrid%perform_fft(forward=.false.)
172 2394138 : mx=real(fftgrid%grid)
173 2394138 : my=aimag(fftgrid%grid) ! TODO: There is NO magic minus here. This is so scuffed....
174 :
175 2394042 : DO i=0,MIN(SIZE(rho,1),size(mx))-1
176 2393946 : rhotot= 0.5*( rho(i,1) + rho(i,2) )
177 2393946 : magmom(i)= SQRT( (0.5*(rho(i,1)-rho(i,2)))**2 + mx(i)**2 + my(i)**2 )
178 3040794 : IF (l_rdm.AND.dograds) rhodiag(i,1:2) = rho(i,1:2)
179 2393946 : rho(i,1)= rhotot+magmom(i)
180 2393946 : rho(i,2)= rhotot-magmom(i)
181 2394042 : IF (l_rdm.AND.dograds) THEN ! prepare rotation matrix
182 323424 : mmx = 2.0 * mx(i)
183 323424 : mmy = 2.0 * my(i)
184 323424 : mmz = rhodiag(i,1) - rhodiag(i,2)
185 323424 : CALL pol_angle(mmx,mmy,mmz,theta,phi)
186 323424 : cossqu(i) = cos(0.5*theta)**2
187 323424 : sinsqu(i) = sin(0.5*theta)**2
188 323424 : sincos(i) = 2.0*sin(0.5*theta)*cos(0.5*theta)
189 323424 : exi(i) = exp(-ImagUnit*phi)
190 : ENDIF
191 : END DO
192 : ENDIF
193 : ENDIF
194 351 : IF (dograds) THEN
195 284 : IF (PRESENT(xcpot)) THEN
196 284 : CALL xcpot%alloc_gradients(griddim,jspins,grad)
197 : END IF
198 :
199 :
200 : ! in non-collinear calculations the derivatives of |m| are calculated in real space.
201 :
202 :
203 284 : IF (.not.l_noco) THEN
204 :
205 : ! In collinear calculations all derivatives are calculated in g-spce,
206 : ndm = 0
207 868 : DO idm = 1,3
208 651 : fd=0.0;fd(idm)=1
209 1572 : DO js=1,jspins
210 921 : if (gmaxz>0.0) then
211 0 : call fftgrid%putFieldOnGrid(stars,den_pw(:,js),cell,gmax,gmaxz,firstderiv=fd)
212 : else
213 921 : call fftgrid%putFieldOnGrid(stars,den_pw(:,js),cell,gmax,firstderiv=fd)
214 : end if
215 921 : call fftgrid%perform_fft(forward=.false.)
216 9179895 : rhd1(0:,js,idm)=fftgrid%grid
217 : END DO
218 868 : IF (allocated(grad%laplace).or.allocated(grad%agrt)) THEN
219 : !Higher derivatives needed
220 1953 : DO jdm = 1,idm
221 1302 : sd=0;sd(jdm)=1
222 1302 : ndm = ndm + 1
223 3795 : DO js=1,jspins
224 1842 : if (gmaxz>0.0) then
225 0 : call fftgrid%putFieldOnGrid(stars,den_pw(:,js),cell,gmax,gmaxz,firstderiv=fd,secondderiv=sd)
226 : else
227 1842 : call fftgrid%putFieldOnGrid(stars,den_pw(:,js),cell,gmax,firstderiv=fd,secondderiv=sd)
228 : end if
229 1842 : call fftgrid%perform_fft(forward=.false.)
230 18359790 : rhd2(0:,js,ndm)=fftgrid%grid
231 : END DO
232 : END DO ! jdm
233 : ENDIF
234 : END DO ! idm
235 :
236 : ELSE !noco case
237 :
238 : IF (l_rdm) THEN
239 :
240 67 : CALL grdrsis( rhodiag(:,1),cell,fftgrid%dimensions,der(:,:,1) )
241 67 : CALL grdrsis( rhodiag(:,2),cell,fftgrid%dimensions,der(:,:,2) )
242 67 : CALL grdrsis( mx,cell,fftgrid%dimensions,der(:,:,3) )
243 67 : CALL grdrsis( my,cell,fftgrid%dimensions,der(:,:,4) )
244 :
245 323491 : DO i=0,griddim-1 ! project on magnetization axis
246 1293763 : DO idm=1,3
247 970272 : rho21 = der(i,idm,3) + ImagUnit * der(i,idm,4)
248 : rhd1(i,1,idm) = cossqu(i) * der(i,idm,1) + &
249 : sincos(i) * real( exi(i)*rho21 ) + &
250 970272 : sinsqu(i) * der(i,idm,2)
251 : rhd1(i,2,idm) = sinsqu(i) * der(i,idm,1) - &
252 : sincos(i) * real( exi(i)*rho21 ) + &
253 1293696 : cossqu(i) * der(i,idm,2)
254 : ENDDO
255 : ENDDO
256 : !Now second derivatives
257 268 : DO idm = 1,3
258 201 : CALL grdrsis( der(:,idm,1),cell,fftgrid%dimensions, dder(:,:,idm,1) )
259 201 : CALL grdrsis( der(:,idm,2),cell,fftgrid%dimensions, dder(:,:,idm,2) )
260 201 : CALL grdrsis( der(:,idm,3),cell,fftgrid%dimensions, dder(:,:,idm,3) )
261 201 : CALL grdrsis( der(:,idm,4),cell,fftgrid%dimensions, dder(:,:,idm,4) )
262 970540 : DO i=0,griddim-1 ! project on magnetization axis
263 3881289 : DO jdm=1,3
264 2910816 : rho21 = dder(i,jdm,idm,3) + ImagUnit * dder(i,jdm,idm,4)
265 : rhdd(i,1,jdm,idm) = cossqu(i) * dder(i,jdm,idm,1) + &
266 : sincos(i) * real( exi(i)*rho21 ) + &
267 2910816 : sinsqu(i) * dder(i,jdm,idm,2)
268 : rhdd(i,2,jdm,idm) = sinsqu(i) * dder(i,jdm,idm,1) - &
269 : sincos(i) * real( exi(i)*rho21 ) + &
270 3881088 : cossqu(i) * dder(i,jdm,idm,2)
271 : ENDDO
272 : ENDDO
273 : ENDDO
274 201 : DO j=1,2
275 647049 : DO i=0,griddim-1
276 646848 : rhd2(i,j,1) = rhdd(i,j,1,1)
277 646848 : rhd2(i,j,2) = 0.5*(rhdd(i,j,1,2)+rhdd(i,j,2,1)) ! xy - averaging should be unneccessary
278 646848 : rhd2(i,j,3) = rhdd(i,j,2,2)
279 646848 : rhd2(i,j,4) = 0.5*(rhdd(i,j,1,3)+rhdd(i,j,3,1)) ! zx
280 646848 : rhd2(i,j,5) = 0.5*(rhdd(i,j,2,3)+rhdd(i,j,3,2)) ! yz
281 646982 : rhd2(i,j,6) = rhdd(i,j,3,3)
282 : ENDDO
283 : ENDDO
284 67 : DEALLOCATE (rhodiag,der,rhdd,dder,sinsqu,cossqu,sincos,exi)
285 :
286 : ELSE
287 :
288 : CALL grdrsis(magmom,cell,fftgrid%dimensions,dmagmom )
289 :
290 : DO i=0,griddim-1
291 : DO idm=1,3
292 : rhotot= rhd1(i,1,idm)/2.+rhd1(i,2,idm)/2.
293 : rhd1(i,1,idm)= rhotot+dmagmom(i,idm)
294 : rhd1(i,2,idm)= rhotot-dmagmom(i,idm)
295 : END DO
296 : END DO
297 : !Second derivatives
298 : DO idm = 1,3
299 : CALL grdrsis(dmagmom(:,idm),cell,fftgrid%dimensions,ddmagmom(:,:,idm) )
300 : END DO
301 : ndm= 0
302 : DO idm = 1,3
303 : DO jdm = 1,idm
304 : ndm = ndm + 1
305 : DO i=0,griddim-1
306 : rhotot= rhd2(i,1,ndm)/2.+rhd2(i,2,ndm)/2.
307 : rhd2(i,1,ndm)= rhotot + ( ddmagmom(i,jdm,idm) + ddmagmom(i,idm,jdm) )/2.
308 : rhd2(i,2,ndm)= rhotot - ( ddmagmom(i,jdm,idm) + ddmagmom(i,idm,jdm) )/2.
309 : END DO
310 : ENDDO !jdm
311 : ENDDO !idm
312 : DEALLOCATE(dmagmom,ddmagmom)
313 : END IF
314 :
315 : END IF
316 :
317 :
318 :
319 :
320 : !
321 : ! calculate the quantities such as abs(grad(rho)),.. used in
322 : ! evaluating the gradient contributions to potential and energy.
323 : !
324 284 : IF (PRESENT(rho)) THEN
325 : CALL mkgxyz3 (rho,rhd1(0:,:,1),rhd1(0:,:,2),rhd1(0:,:,3),&
326 276 : rhd2(0:,:,1),rhd2(0:,:,3),rhd2(0:,:,6), rhd2(0:,:,5),rhd2(0:,:,4),rhd2(0:,:,2),0,grad)
327 : ELSE
328 : !Dummy rho (only possible if grad is used for libxc mode)
329 : !CALL mkgxyz3 (RESHAPE((/0.0/),(/1,1/)),rhd1(0:,:,1),rhd1(0:,:,2),rhd1(0:,:,3),&
330 : ! rhd2(0:,:,1),rhd2(0:,:,3),rhd2(0:,:,6), rhd2(0:,:,5),rhd2(0:,:,4),rhd2(0:,:,2),grad)
331 :
332 8 : IF (dograds.and.(.not.PRESENT(xcpot))) THEN
333 0 : ALLOCATE(grad%gr(3,griddim,1))
334 : END IF
335 :
336 : CALL mkgxyz3 (0*rhd1(0:,:,1),rhd1(0:,:,1),rhd1(0:,:,2),rhd1(0:,:,3),&
337 98236 : rhd2(0:,:,1),rhd2(0:,:,3),rhd2(0:,:,6), rhd2(0:,:,5),rhd2(0:,:,4),rhd2(0:,:,2),0,grad)
338 : END IF
339 :
340 : ENDIF
341 351 : rhoCutGrad = 3.0e-6
342 351 : IF (PRESENT(rho)) THEN
343 10086206 : WHERE(ABS(rho) < d_15) rho = d_15
344 343 : IF(PRESENT(xcpot)) THEN
345 : SELECT TYPE(xcpot)
346 : TYPE IS (t_xcpot_inbuild)
347 332 : IF (dograds) THEN
348 : ! In the following, gradients at places with very small densities are set to zero. This is done for stability reasons.
349 2770113 : DO i = 1, griddim
350 2769841 : IF(rho(i-1,1).LT.rhoCutGrad) THEN
351 7060 : grad%agru(i) = 0.0
352 7060 : grad%g2ru(i) = 0.0
353 7060 : grad%gggru(i) = 0.0
354 :
355 7060 : grad%agrt(i) = grad%agrd(i)
356 7060 : grad%g2rt(i) = grad%g2rd(i)
357 7060 : grad%gggrt(i) = grad%gggrd(i)
358 : END IF
359 :
360 2769841 : IF(rho(i-1,jspins).LT.rhoCutGrad) THEN
361 7256 : grad%agrd(i) = 0.0
362 7256 : grad%g2rd(i) = 0.0
363 7256 : grad%gggrd(i) = 0.0
364 :
365 7256 : grad%agrt(i) = grad%agru(i)
366 7256 : grad%g2rt(i) = grad%g2ru(i)
367 7256 : grad%gggrt(i) = grad%gggru(i)
368 : END IF
369 :
370 2770113 : IF((rho(i-1,1).LT.rhoCutGrad).AND.(rho(i-1,jspins).LT.rhoCutGrad)) THEN
371 7060 : grad%gzgr(i) = 0.0
372 : END IF
373 : END DO
374 : END IF
375 : END SELECT
376 : END IF
377 : ENDIF
378 351 : IF (PRESENT(rhoim)) THEN
379 0 : WHERE(ABS(rhoim) < d_15) rhoim = d_15
380 : ENDIF
381 :
382 351 : END SUBROUTINE pw_to_grid
383 :
384 :
385 1362 : SUBROUTINE pw_from_grid(stars,v_in,v_out_pw,v_out_pw_w)
386 : USE m_convol
387 : USE m_types
388 : IMPLICIT NONE
389 : TYPE(t_stars),INTENT(IN) :: stars
390 : REAL,INTENT(INOUT) :: v_in(0:,:)
391 : COMPLEX,INTENT(INOUT) :: v_out_pw(:,:)
392 : COMPLEX,INTENT(INOUT),OPTIONAL:: v_out_pw_w(:,:)
393 :
394 :
395 : INTEGER :: js,k,i
396 1362 : REAL,ALLOCATABLE :: bf3(:),vcon(:)
397 1362 : COMPLEX, ALLOCATABLE :: fg3(:)
398 1362 : if (present(v_out_pw_w)) ALLOCATE( bf3(size(stars%ufft)),vcon(size(stars%ufft)))
399 1362 : ALLOCATE ( fg3(stars%ng3) )
400 3333 : DO js = 1,SIZE(v_in,2)
401 36536920 : fftgrid%grid=v_in(0:,js)
402 1971 : call fftgrid%perform_fft(forward=.true.)
403 1971 : if (gmaxz>0.0) then
404 0 : call fftgrid%takeFieldFromGrid(stars,fg3,gmax,gmaxz)
405 : else
406 1971 : call fftgrid%takeFieldFromGrid(stars,fg3,gmax)
407 : END IF
408 4346987 : v_out_pw(:,js) = v_out_pw(:,js) + fg3(:)
409 :
410 : !----> add to warped coulomb potential
411 3333 : IF (present(v_out_pw_w)) THEN
412 1415 : if (size(fftgrid%grid)==size(stars%ufft)) THEN
413 16338444 : fftgrid%grid=v_in(0:,js)*stars%ufft
414 291 : call fftgrid%perform_fft(forward=.true.)
415 291 : if (gmaxz>0.0) then
416 0 : call fftgrid%takeFieldFromGrid(stars,fg3,gmax,gmaxz)
417 : else
418 291 : call fftgrid%takeFieldFromGrid(stars,fg3,gmax)
419 : end if
420 1545476 : fg3 = fg3*stars%nstr
421 : else
422 1124 : call convol(stars,fg3)
423 : ENDIF
424 3148886 : v_out_pw_w(:,js) = v_out_pw_w(:,js) + fg3
425 : ENDIF
426 : END DO
427 1362 : END SUBROUTINE pw_from_grid
428 :
429 339 : SUBROUTINE finish_pw_grid()
430 : IMPLICIT NONE
431 : !IF (ALLOCATED(igxc_fft)) DEALLOCATE(igxc_fft,gxc_fft)
432 339 : END SUBROUTINE finish_pw_grid
433 :
434 337 : END MODULE m_pw_tofrom_grid
|
