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 590 : 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 590 : if (present(xcpot)) THEN
36 588 : gmax=xcpot%gmaxxc
37 588 : gmaxz=stars%gmaxz
38 588 : if (xcpot%needs_grad().and.gmax>0.0) then
39 258 : if (gmaxz>0.0) then
40 0 : call fftgrid%init(cell,sym,gmax,gmaxz)
41 : else
42 258 : call fftgrid%init(cell,sym,gmax)
43 : end if
44 : else
45 1320 : call fftgrid%init((/3*stars%mx1,3*stars%mx2,3*stars%mx3/))
46 330 : 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 590 : griddim=size(fftgrid%grid)
58 :
59 590 : END SUBROUTINE init_pw_grid
60 :
61 1682 : 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 841 : COMPLEX, ALLOCATABLE :: cqpw(:,:),ph_wrk(:)
114 841 : REAL, ALLOCATABLE :: bf3(:)
115 841 : REAL, ALLOCATABLE :: rhd1(:,:,:),rhd2(:,:,:)
116 841 : REAL, ALLOCATABLE :: mx(:),my(:)
117 841 : REAL, ALLOCATABLE :: magmom(:),dmagmom(:,:),ddmagmom(:,:,:)
118 841 : REAL, ALLOCATABLE :: rhodiag(:,:),der(:,:,:),dder(:,:,:,:)
119 841 : REAL, ALLOCATABLE :: sinsqu(:),cossqu(:),sincos(:),rhdd(:,:,:,:)
120 841 : COMPLEX, ALLOCATABLE :: exi(:)
121 :
122 : LOGICAL, PARAMETER :: l_rdm=.true.
123 :
124 : ! Allocate arrays
125 841 : ALLOCATE( bf3(0:griddim-1))
126 841 : IF (dograds) THEN
127 524 : IF (PRESENT(rho)) ALLOCATE(rho(0:griddim-1,jspins))
128 266 : IF (PRESENT(rhoim)) ALLOCATE(rhoim(0:griddim-1,jspins))
129 798 : ALLOCATE( ph_wrk(0:griddim-1),rhd1(0:griddim-1,jspins,3))
130 798 : ALLOCATE( rhd2(0:griddim-1,jspins,6) )
131 : ELSE
132 1150 : IF (PRESENT(rho)) ALLOCATE(rho(0:griddim-1,jspins))
133 814 : IF (PRESENT(rhoim)) ALLOCATE(rhoim(0:griddim-1,jspins))
134 : ENDIF
135 841 : IF (l_noco) THEN
136 92 : IF (dograds) THEN
137 65 : ALLOCATE( mx(0:griddim-1),my(0:griddim-1),magmom(0:griddim-1))
138 : IF (l_rdm) THEN
139 130 : 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 65 : 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 27 : 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 841 : IF (PRESENT(rho)) THEN
150 : !Put den_pw on grid and store into rho(:,1:2)
151 1844 : DO js=1,jspins
152 1011 : if (gmaxz>0.0) then
153 0 : call fftgrid%putFieldOnGrid(stars,den_pw(:,js),cell,gmax,gmaxz)
154 : else
155 1011 : call fftgrid%putFieldOnGrid(stars,den_pw(:,js),cell,gmax)
156 : end if
157 1011 : 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 14872971 : rho(0:,js) = REAL(fftgrid%grid)
161 2506094 : IF (PRESENT(rhoim)) rhoim(0:,js) = AIMAG(fftgrid%grid)
162 : END DO
163 :
164 833 : IF (l_noco) THEN
165 : ! Get mx,my on real space grid and recalculate rho and magmom
166 92 : if (gmaxz>0.0) then
167 0 : call fftgrid%putFieldOnGrid(stars,den_pw(:,3),cell,gmax,gmaxz)
168 : else
169 92 : call fftgrid%putFieldOnGrid(stars,den_pw(:,3),cell,gmax)
170 : end if
171 92 : call fftgrid%perform_fft(forward=.false.)
172 2289228 : mx=real(fftgrid%grid)
173 2289228 : my=aimag(fftgrid%grid) ! TODO: There is NO magic minus here. This is so scuffed....
174 :
175 2289136 : DO i=0,MIN(SIZE(rho,1),size(mx))-1
176 2289044 : rhotot= 0.5*( rho(i,1) + rho(i,2) )
177 2289044 : magmom(i)= SQRT( (0.5*(rho(i,1)-rho(i,2)))**2 + mx(i)**2 + my(i)**2 )
178 2804820 : IF (l_rdm.AND.dograds) rhodiag(i,1:2) = rho(i,1:2)
179 2289044 : rho(i,1)= rhotot+magmom(i)
180 2289044 : rho(i,2)= rhotot-magmom(i)
181 2289136 : IF (l_rdm.AND.dograds) THEN ! prepare rotation matrix
182 257888 : mmx = 2.0 * mx(i)
183 257888 : mmy = 2.0 * my(i)
184 257888 : mmz = rhodiag(i,1) - rhodiag(i,2)
185 257888 : CALL pol_angle(mmx,mmy,mmz,theta,phi)
186 257888 : cossqu(i) = cos(0.5*theta)**2
187 257888 : sinsqu(i) = sin(0.5*theta)**2
188 257888 : sincos(i) = 2.0*sin(0.5*theta)*cos(0.5*theta)
189 257888 : exi(i) = exp(-ImagUnit*phi)
190 : ENDIF
191 : END DO
192 : ENDIF
193 : ENDIF
194 841 : IF (dograds) THEN
195 266 : IF (PRESENT(xcpot)) THEN
196 266 : 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 266 : IF (.not.l_noco) THEN
204 :
205 : ! In collinear calculations all derivatives are calculated in g-spce,
206 : ndm = 0
207 804 : DO idm = 1,3
208 603 : fd=0.0;fd(idm)=1
209 1428 : DO js=1,jspins
210 825 : if (gmaxz>0.0) then
211 0 : call fftgrid%putFieldOnGrid(stars,den_pw(:,js),cell,gmax,gmaxz,firstderiv=fd)
212 : else
213 825 : call fftgrid%putFieldOnGrid(stars,den_pw(:,js),cell,gmax,firstderiv=fd)
214 : end if
215 825 : call fftgrid%perform_fft(forward=.false.)
216 8074947 : rhd1(0:,js,idm)=fftgrid%grid
217 : END DO
218 804 : IF (allocated(grad%laplace).or.allocated(grad%agrt)) THEN
219 : !Higher derivatives needed
220 1809 : DO jdm = 1,idm
221 1206 : sd=0;sd(jdm)=1
222 1206 : ndm = ndm + 1
223 3459 : DO js=1,jspins
224 1650 : if (gmaxz>0.0) then
225 0 : call fftgrid%putFieldOnGrid(stars,den_pw(:,js),cell,gmax,gmaxz,firstderiv=fd,secondderiv=sd)
226 : else
227 1650 : call fftgrid%putFieldOnGrid(stars,den_pw(:,js),cell,gmax,firstderiv=fd,secondderiv=sd)
228 : end if
229 1650 : call fftgrid%perform_fft(forward=.false.)
230 16149894 : 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 65 : CALL grdrsis( rhodiag(:,1),cell,fftgrid%dimensions,der(:,:,1) )
241 65 : CALL grdrsis( rhodiag(:,2),cell,fftgrid%dimensions,der(:,:,2) )
242 65 : CALL grdrsis( mx,cell,fftgrid%dimensions,der(:,:,3) )
243 65 : CALL grdrsis( my,cell,fftgrid%dimensions,der(:,:,4) )
244 :
245 257953 : DO i=0,griddim-1 ! project on magnetization axis
246 1031617 : DO idm=1,3
247 773664 : 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 773664 : 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 1031552 : cossqu(i) * der(i,idm,2)
254 : ENDDO
255 : ENDDO
256 : !Now second derivatives
257 260 : DO idm = 1,3
258 195 : CALL grdrsis( der(:,idm,1),cell,fftgrid%dimensions, dder(:,:,idm,1) )
259 195 : CALL grdrsis( der(:,idm,2),cell,fftgrid%dimensions, dder(:,:,idm,2) )
260 195 : CALL grdrsis( der(:,idm,3),cell,fftgrid%dimensions, dder(:,:,idm,3) )
261 195 : CALL grdrsis( der(:,idm,4),cell,fftgrid%dimensions, dder(:,:,idm,4) )
262 773924 : DO i=0,griddim-1 ! project on magnetization axis
263 3094851 : DO jdm=1,3
264 2320992 : 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 2320992 : 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 3094656 : cossqu(i) * dder(i,jdm,idm,2)
271 : ENDDO
272 : ENDDO
273 : ENDDO
274 195 : DO j=1,2
275 515971 : DO i=0,griddim-1
276 515776 : rhd2(i,j,1) = rhdd(i,j,1,1)
277 515776 : rhd2(i,j,2) = 0.5*(rhdd(i,j,1,2)+rhdd(i,j,2,1)) ! xy - averaging should be unneccessary
278 515776 : rhd2(i,j,3) = rhdd(i,j,2,2)
279 515776 : rhd2(i,j,4) = 0.5*(rhdd(i,j,1,3)+rhdd(i,j,3,1)) ! zx
280 515776 : rhd2(i,j,5) = 0.5*(rhdd(i,j,2,3)+rhdd(i,j,3,2)) ! yz
281 515906 : rhd2(i,j,6) = rhdd(i,j,3,3)
282 : ENDDO
283 : ENDDO
284 65 : 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 266 : IF (PRESENT(rho)) THEN
325 : CALL mkgxyz3 (rho,rhd1(0:,:,1),rhd1(0:,:,2),rhd1(0:,:,3),&
326 258 : 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 841 : rhoCutGrad = 3.0e-6
342 841 : IF (PRESENT(rho)) THEN
343 14873804 : WHERE(ABS(rho) < d_15) rho = d_15
344 833 : IF(PRESENT(xcpot)) THEN
345 : SELECT TYPE(xcpot)
346 : TYPE IS (t_xcpot_inbuild)
347 312 : 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 2520425 : DO i = 1, griddim
350 2520171 : IF(rho(i-1,1).LT.rhoCutGrad) THEN
351 7153 : grad%agru(i) = 0.0
352 7153 : grad%g2ru(i) = 0.0
353 7153 : grad%gggru(i) = 0.0
354 :
355 7153 : grad%agrt(i) = grad%agrd(i)
356 7153 : grad%g2rt(i) = grad%g2rd(i)
357 7153 : grad%gggrt(i) = grad%gggrd(i)
358 : END IF
359 :
360 2520171 : IF(rho(i-1,jspins).LT.rhoCutGrad) THEN
361 7487 : grad%agrd(i) = 0.0
362 7487 : grad%g2rd(i) = 0.0
363 7487 : grad%gggrd(i) = 0.0
364 :
365 7487 : grad%agrt(i) = grad%agru(i)
366 7487 : grad%g2rt(i) = grad%g2ru(i)
367 7487 : grad%gggrt(i) = grad%gggru(i)
368 : END IF
369 :
370 2520425 : IF((rho(i-1,1).LT.rhoCutGrad).AND.(rho(i-1,jspins).LT.rhoCutGrad)) THEN
371 7153 : grad%gzgr(i) = 0.0
372 : END IF
373 : END DO
374 : END IF
375 : END SELECT
376 : END IF
377 : ENDIF
378 841 : IF (PRESENT(rhoim)) THEN
379 2504728 : WHERE(ABS(rhoim) < d_15) rhoim = d_15
380 : ENDIF
381 :
382 841 : END SUBROUTINE pw_to_grid
383 :
384 :
385 1888 : 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 1888 : REAL,ALLOCATABLE :: bf3(:),vcon(:)
397 1888 : COMPLEX, ALLOCATABLE :: fg3(:)
398 1888 : if (present(v_out_pw_w)) ALLOCATE( bf3(size(stars%ufft)),vcon(size(stars%ufft)))
399 1888 : ALLOCATE ( fg3(stars%ng3) )
400 4322 : DO js = 1,SIZE(v_in,2)
401 40974938 : fftgrid%grid=v_in(0:,js)
402 2434 : call fftgrid%perform_fft(forward=.true.)
403 2434 : if (gmaxz>0.0) then
404 0 : call fftgrid%takeFieldFromGrid(stars,fg3,gmax,gmaxz)
405 : else
406 2434 : call fftgrid%takeFieldFromGrid(stars,fg3,gmax)
407 : END IF
408 4545340 : v_out_pw(:,js) = v_out_pw(:,js) + fg3(:)
409 :
410 : !----> add to warped coulomb potential
411 4322 : IF (present(v_out_pw_w)) THEN
412 1409 : if (size(fftgrid%grid)==size(stars%ufft)) THEN
413 17237415 : fftgrid%grid=v_in(0:,js)*stars%ufft
414 375 : call fftgrid%perform_fft(forward=.true.)
415 375 : if (gmaxz>0.0) then
416 0 : call fftgrid%takeFieldFromGrid(stars,fg3,gmax,gmaxz)
417 : else
418 375 : call fftgrid%takeFieldFromGrid(stars,fg3,gmax)
419 : end if
420 1631789 : fg3 = fg3*stars%nstr
421 : else
422 1034 : call convol(stars,fg3)
423 : ENDIF
424 2945645 : v_out_pw_w(:,js) = v_out_pw_w(:,js) + fg3
425 : ENDIF
426 : END DO
427 1888 : END SUBROUTINE pw_from_grid
428 :
429 351 : SUBROUTINE finish_pw_grid()
430 : IMPLICIT NONE
431 : !IF (ALLOCATED(igxc_fft)) DEALLOCATE(igxc_fft,gxc_fft)
432 351 : END SUBROUTINE finish_pw_grid
433 :
434 827 : END MODULE m_pw_tofrom_grid
|
