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 : PRIVATE
9 : REAL,PARAMETER:: d_15=1.e-15
10 :
11 : INTEGER :: ifftd,ifftxc3
12 : !-----> fft information for xc potential + energy
13 : INTEGER, ALLOCATABLE :: igxc_fft(:)
14 : REAL, ALLOCATABLE :: gxc_fft(:,:) !gxc_fft(ig,idm)
15 :
16 : PUBLIC :: init_pw_grid,pw_to_grid,pw_from_grid,finish_pw_grid
17 : CONTAINS
18 170 : SUBROUTINE init_pw_grid(xcpot,stars,sym,cell)
19 : USE m_prpxcfftmap
20 : USE m_types
21 : IMPLICIT NONE
22 : CLASS(t_xcpot),INTENT(IN) :: xcpot
23 : TYPE(t_stars),INTENT(IN) :: stars
24 : TYPE(t_sym),INTENT(IN) :: sym
25 : TYPE(t_cell),INTENT(IN) :: cell
26 :
27 : !---> set up pointer for backtransformation of from g-vector in
28 : ! positive domain of xc density fftbox into stars.
29 : ! also the x,y,z components of the g-vectors are set up to calculate
30 : ! derivatives.
31 : ! in principle this can also be done in main program once.
32 : ! it is done here to save memory.
33 : !
34 :
35 170 : ifftd=27*stars%mx1*stars%mx2*stars%mx3
36 170 : ifftxc3 = stars%kxc1_fft*stars%kxc2_fft*stars%kxc3_fft
37 170 : IF (xcpot%needs_grad()) THEN
38 161 : CALL prp_xcfft_map(stars,sym, cell, igxc_fft,gxc_fft)
39 : ENDIF
40 :
41 170 : END SUBROUTINE init_pw_grid
42 :
43 380 : SUBROUTINE pw_to_grid(xcpot,jspins,l_noco,stars,cell,den_pw,grad,rho)
44 : !.....------------------------------------------------------------------
45 : !-------> abbreviations
46 : !
47 : ! ph_wrk: work array containing phase * g_x,gy......
48 : ! den%pw: charge density stored as stars
49 : ! rho : charge density stored in real space
50 : ! v_xc : exchange-correlation potential in real space
51 : ! exc : exchange-correlation energy density in real space
52 : ! kxc1d : dimension of the charge density fft box in the pos. domain
53 : ! kxc2d : defined in dimens.f program (subroutine apws).1,2,3 indic
54 : ! kxc3d ; a_1, a_2, a_3 directions.
55 : ! kq(i) : i=1,2,3 actual length of the fft-box for which fft is done
56 : ! nstr : number of members (arms) of reciprocal lattice (g) vector
57 : ! of each star.
58 : ! nxc3_fft: number of stars in the charge density fft-box
59 : ! ng3 : number of 3 dim. stars in the charge density sphere define
60 : ! by gmax
61 : ! kmxxc_fft: number of g-vectors forming the nxc3_fft stars in the
62 : ! charge density or xc-density sphere
63 : ! kimax : number of g-vectors forming the ng3 stars in the gmax-sphe
64 : ! igfft : pointer from the g-sphere (stored as stars) to fft-grid
65 : ! and from fft-grid to g-sphere (stored as stars)
66 : ! pgfft : contains the phases of the g-vectors of sph.
67 : ! isn : isn = +1, fft transform for g-space to r-space
68 : ! isn = -1, vice versa
69 : !
70 : !-------------------------------------------------------------------
71 : USE m_grdrsis
72 : USE m_mkgxyz3
73 : USE m_fft3dxc
74 : USE m_fft3d
75 : USE m_types
76 : use m_constants
77 : IMPLICIT NONE
78 : CLASS(t_xcpot),INTENT(IN) :: xcpot
79 : INTEGER,INTENT(IN) :: jspins
80 : LOGICAL,INTENT(IN) :: l_noco
81 : TYPE(t_stars),INTENT(IN) :: stars
82 : TYPE(t_cell),INTENT(IN) :: cell
83 : COMPLEX,INTENT(IN) :: den_pw(:,:)
84 : TYPE(t_gradients),INTENT(OUT) :: grad
85 : REAL,ALLOCATABLE,INTENT(out),OPTIONAL :: rho(:,:)
86 :
87 :
88 : INTEGER :: js,i,idm,ig,ndm,jdm
89 : REAL :: rhotot
90 : ! .. Local Arrays ..
91 380 : COMPLEX, ALLOCATABLE :: cqpw(:,:),ph_wrk(:)
92 190 : REAL, ALLOCATABLE :: bf3(:)
93 190 : REAL, ALLOCATABLE :: rhd1(:,:,:),rhd2(:,:,:)
94 190 : REAL, ALLOCATABLE :: mx(:),my(:)
95 190 : REAL, ALLOCATABLE :: magmom(:),dmagmom(:,:),ddmagmom(:,:,:)
96 :
97 : ! Allocate arrays
98 190 : ALLOCATE( bf3(0:ifftd-1))
99 190 : IF (xcpot%needs_grad()) THEN
100 181 : IF (PRESENT(rho)) ALLOCATE(rho(0:ifftxc3-1,jspins))
101 181 : ALLOCATE( ph_wrk(0:ifftxc3-1),rhd1(0:ifftxc3-1,jspins,3))
102 181 : ALLOCATE( rhd2(0:ifftxc3-1,jspins,6) )
103 : ELSE
104 9 : IF (PRESENT(rho)) ALLOCATE(rho(0:ifftd-1,jspins))
105 : ENDIF
106 190 : IF (l_noco) THEN
107 124 : IF (xcpot%needs_grad()) THEN
108 124 : ALLOCATE( mx(0:ifftxc3-1),my(0:ifftxc3-1),magmom(0:ifftxc3-1))
109 124 : ALLOCATE(dmagmom(0:ifftxc3-1,3),ddmagmom(0:ifftxc3-1,3,3) )
110 : ELSE
111 0 : ALLOCATE( mx(0:ifftd-1),my(0:ifftd-1),magmom(0:ifftd-1))
112 : ENDIF
113 : END IF
114 :
115 190 : IF (PRESENT(rho)) THEN
116 : !Put den_pw on grid and store into rho(:,1:2)
117 474 : DO js=1,jspins
118 474 : IF (xcpot%needs_grad()) THEN
119 : CALL fft3dxc(rho(0:,js),bf3, den_pw(:,js), stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,&
120 295 : stars%nxc3_fft,stars%kmxxc_fft,+1, stars%igfft(0:,1),igxc_fft,stars%pgfft,stars%nstr)
121 : ELSE
122 9 : CALL fft3d(rho(0,js),bf3, den_pw(:,js), stars,+1)
123 : ENDIF
124 : END DO
125 :
126 170 : IF (l_noco) THEN
127 : ! Get mx,my on real space grid and recalculate rho and magmom
128 124 : IF (xcpot%needs_grad()) THEN
129 : CALL fft3dxc(mx,my, den_pw(:,3), stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,&
130 124 : stars%nxc3_fft,stars%kmxxc_fft,+1, stars%igfft(0:,1),igxc_fft,stars%pgfft,stars%nstr)
131 : ELSE
132 0 : CALL fft3d(mx,my, den_pw(:,3), stars,+1)
133 : ENDIF
134 614524 : DO i=0,MIN(SIZE(rho,1),size(mx))-1
135 614400 : rhotot= 0.5*( rho(i,1) + rho(i,2) )
136 614400 : magmom(i)= SQRT( (0.5*(rho(i,1)-rho(i,2)))**2 + mx(i)**2 + my(i)**2 )
137 614400 : rho(i,1)= rhotot+magmom(i)
138 614524 : rho(i,2)= rhotot-magmom(i)
139 : END DO
140 : ENDIF
141 : ENDIF
142 190 : IF (xcpot%needs_grad()) THEN
143 :
144 : ! In collinear calculations all derivatives are calculated in g-spce,
145 : ! in non-collinear calculations the derivatives of |m| are calculated in real space.
146 :
147 : !--> for d(rho)/d(x,y,z) = rhd1(:,:,idm) (idm=1,2,3).
148 : !
149 : ! ph_wrk: exp(i*(g_x,g_y,g_z)*tau) * g_(x,y,z).
150 :
151 181 : ALLOCATE(cqpw(stars%ng3,jspins))
152 :
153 536 : cqpw(:,:)= ImagUnit*den_pw(:,:jspins)
154 :
155 1267 : DO idm=1,3
156 1313202 : DO ig = 0 , stars%kmxxc_fft - 1
157 1313202 : ph_wrk(ig) = stars%pgfft(ig) * gxc_fft(ig,idm)
158 : END DO
159 :
160 2854 : DO js=1,jspins
161 : CALL fft3dxc(rhd1(0:,js,idm),bf3, cqpw(:,js), stars%kxc1_fft,stars%kxc2_fft,&
162 1608 : stars%kxc3_fft,stars%nxc3_fft,stars%kmxxc_fft,+1, stars%igfft(0:,1),igxc_fft,ph_wrk,stars%nstr)
163 : END DO
164 : END DO
165 :
166 181 : IF (l_noco) THEN
167 :
168 124 : CALL grdrsis(magmom,cell,stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,dmagmom )
169 :
170 614524 : DO i=0,ifftxc3-1
171 1843324 : DO idm=1,3
172 1843200 : rhotot= rhd1(i,1,idm)/2.+rhd1(i,2,idm)/2.
173 1843200 : rhd1(i,1,idm)= rhotot+dmagmom(i,idm)
174 2457600 : rhd1(i,2,idm)= rhotot-dmagmom(i,idm)
175 : END DO
176 : END DO
177 : END IF
178 :
179 : !--> for dd(rho)/d(xx,xy,yy,zx,yz,zz) = rhd2(:,:,idm) (idm=1,2,3,4,5,6)
180 : !
181 : ! ph_wrk: exp(i*(g_x,g_y,g_z)*tau) * g_(x,y,z) * g_(x,y,z)
182 :
183 536 : cqpw(:,:)= -den_pw(:,:jspins)
184 :
185 : ndm = 0
186 1267 : DO idm = 1,3
187 2896 : DO jdm = 1,idm
188 1086 : ndm = ndm + 1
189 2626404 : DO ig = 0 , stars%kmxxc_fft-1
190 2626404 : ph_wrk(ig) = stars%pgfft(ig)*gxc_fft(ig,idm)*gxc_fft(ig,jdm)
191 : ENDDO
192 :
193 5889 : DO js=1,jspins
194 : CALL fft3dxc(rhd2(0:,js,ndm),bf3, cqpw(:,js), stars%kxc1_fft,stars%kxc2_fft,&
195 3216 : stars%kxc3_fft,stars%nxc3_fft,stars%kmxxc_fft,+1, stars%igfft(0:,1),igxc_fft,ph_wrk,stars%nstr)
196 : END DO
197 : END DO ! jdm
198 : END DO ! idm
199 :
200 181 : DEALLOCATE(cqpw)
201 :
202 181 : IF (l_noco) THEN
203 868 : DO idm = 1,3
204 496 : CALL grdrsis(dmagmom(0,idm),cell,stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,ddmagmom(0,1,idm) )
205 : END DO
206 : ndm= 0
207 868 : DO idm = 1,3
208 1984 : DO jdm = 1,idm
209 744 : ndm = ndm + 1
210 3687516 : DO i=0,ifftxc3-1
211 3686400 : rhotot= rhd2(i,1,ndm)/2.+rhd2(i,2,ndm)/2.
212 3686400 : rhd2(i,1,ndm)= rhotot + ( ddmagmom(i,jdm,idm) + ddmagmom(i,idm,jdm) )/2.
213 3687144 : rhd2(i,2,ndm)= rhotot - ( ddmagmom(i,jdm,idm) + ddmagmom(i,idm,jdm) )/2.
214 : END DO
215 : ENDDO !jdm
216 : ENDDO !idm
217 : END IF
218 181 : CALL xcpot%alloc_gradients(ifftxc3,jspins,grad)
219 :
220 : !
221 : ! calculate the quantities such as abs(grad(rho)),.. used in
222 : ! evaluating the gradient contributions to potential and energy.
223 : !
224 181 : IF (PRESENT(rho)) THEN
225 : CALL mkgxyz3 (rho,rhd1(0:,:,1),rhd1(0:,:,2),rhd1(0:,:,3),&
226 161 : rhd2(0:,:,1),rhd2(0:,:,3),rhd2(0:,:,6), rhd2(0:,:,5),rhd2(0:,:,4),rhd2(0:,:,2),grad)
227 : ELSE
228 : !Dummy rho (only possible if grad is used for libxc mode)
229 : CALL mkgxyz3 (RESHAPE((/0.0/),(/1,1/)),rhd1(0:,:,1),rhd1(0:,:,2),rhd1(0:,:,3),&
230 20 : rhd2(0:,:,1),rhd2(0:,:,3),rhd2(0:,:,6), rhd2(0:,:,5),rhd2(0:,:,4),rhd2(0:,:,2),grad)
231 : END IF
232 :
233 : ENDIF
234 190 : IF (PRESENT(rho)) THEN
235 170 : WHERE(ABS(rho) < d_15) rho = d_15
236 : ENDIF
237 :
238 190 : END SUBROUTINE pw_to_grid
239 :
240 :
241 530 : SUBROUTINE pw_from_grid(xcpot,stars,l_pw_w,v_in,v_out_pw,v_out_pw_w)
242 : USE m_fft3d
243 : USE m_fft3dxc
244 : USE m_types
245 : IMPLICIT NONE
246 : CLASS(t_xcpot),INTENT(in) :: xcpot
247 : TYPE(t_stars),INTENT(IN) :: stars
248 : REAL,INTENT(INOUT) :: v_in(0:,:)
249 : LOGICAL,INTENT(in) :: l_pw_w
250 : COMPLEX,INTENT(INOUT) :: v_out_pw(:,:)
251 : COMPLEX,INTENT(INOUT),OPTIONAL:: v_out_pw_w(:,:)
252 :
253 :
254 : INTEGER :: js,k,i
255 1060 : REAL,ALLOCATABLE :: bf3(:),vcon(:)
256 530 : COMPLEX, ALLOCATABLE :: fg3(:)
257 530 : ALLOCATE( bf3(0:ifftd-1),fg3(stars%ng3))
258 530 : ALLOCATE ( vcon(0:ifftd-1) )
259 1368 : DO js = 1,SIZE(v_in,2)
260 19898569 : bf3=0.0
261 838 : IF (xcpot%needs_grad()) THEN
262 : CALL fft3dxc(v_in(0:,js),bf3, fg3, stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,&
263 811 : stars%nxc3_fft,stars%kmxxc_fft,-1, stars%igfft(0:,1),igxc_fft,stars%pgfft,stars%nstr)
264 : ELSE
265 27 : vcon(0:)=v_in(0:,js)
266 27 : CALL fft3d(v_in(0:,js),bf3, fg3, stars,-1)
267 : ENDIF
268 482439 : DO k = 1,MERGE(stars%nxc3_fft,stars%ng3,xcpot%needs_grad())
269 482439 : v_out_pw(k,js) = v_out_pw(k,js) + fg3(k)
270 : ENDDO
271 :
272 1368 : IF (l_pw_w) THEN
273 778 : IF (xcpot%needs_grad()) THEN
274 : !----> Perform fft transform: v_xc(star) --> vxc(r)
275 : ! !Use large fft mesh for convolution
276 751 : fg3(stars%nxc3_fft+1:)=0.0
277 751 : CALL fft3d(vcon(0),bf3, fg3, stars,+1)
278 : ENDIF
279 : !
280 : !----> Convolute with step function
281 : !
282 19025329 : DO i=0,ifftd-1
283 778 : vcon(i)=stars%ufft(i)*vcon(i)
284 : ENDDO
285 19025329 : bf3=0.0
286 778 : CALL fft3d(vcon(0),bf3, fg3, stars,-1)
287 778 : fg3=fg3*stars%nstr
288 : !
289 : !----> add to warped coulomb potential
290 : !
291 776523 : DO k = 1,stars%ng3
292 776523 : v_out_pw_w(k,js) = v_out_pw_w(k,js) + fg3(k)
293 : ENDDO
294 : ENDIF
295 : END DO
296 530 : END SUBROUTINE pw_from_grid
297 :
298 170 : SUBROUTINE finish_pw_grid()
299 : IMPLICIT NONE
300 170 : IF (ALLOCATED(igxc_fft)) DEALLOCATE(igxc_fft,gxc_fft)
301 170 : END SUBROUTINE finish_pw_grid
302 :
303 : END MODULE m_pw_tofrom_grid
|
