Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2019 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_divergence
7 : USE m_types
8 : USE m_juDFT
9 : PRIVATE
10 : PUBLIC :: divergence, vac_grad, divpotgrad
11 :
12 : CONTAINS
13 64 : SUBROUTINE divergence(input,stars,atoms,sphhar,vacuum,sym,cell,noco,bxc,div)
14 : USE m_lattHarmsSphHarmsConv
15 : USE m_gradYlm
16 : USE m_constants
17 :
18 : !--------------------------------------------------------------------------
19 : ! Use the interstitial/vacuum divergence subroutine and an external MT-gra-
20 : ! dient routine from juPhon to assemble the divergence of a field into a
21 : ! t_potden variable. The MT-gradient is first calculated in sperical har-
22 : ! monics coefficients.
23 : !--------------------------------------------------------------------------
24 :
25 : IMPLICIT NONE
26 :
27 : TYPE(t_input), INTENT(IN) :: input
28 : TYPE(t_stars), INTENT(IN) :: stars
29 : TYPE(t_atoms), INTENT(IN) :: atoms
30 : TYPE(t_sphhar), INTENT(IN) :: sphhar
31 : TYPE(t_vacuum), INTENT(IN) :: vacuum
32 : TYPE(t_sym), INTENT(IN) :: sym
33 : TYPE(t_cell), INTENT(IN) :: cell
34 : TYPE(t_noco), INTENT(IN) :: noco
35 : TYPE(t_potden), DIMENSION(3), INTENT(INOUT) :: bxc
36 : TYPE(t_potden), INTENT(INOUT) :: div
37 :
38 3392 : TYPE(t_potden), DIMENSION(3) :: grad
39 :
40 : INTEGER :: i,iType,indmax, lh
41 64 : COMPLEX, ALLOCATABLE :: flm(:,:,:),grsflm1(:,:,:,:),grsflm2(:,:,:,:),grsflm3(:,:,:,:),divflm(:,:,:) ! (iR,lm,n[,x,i])
42 :
43 64 : CALL timestart("MT divergence")
44 64 : indmax=(atoms%lmaxd+1)**2
45 :
46 320 : ALLOCATE(flm(atoms%jmtd,indmax,atoms%ntype))
47 256 : ALLOCATE(divflm(atoms%jmtd,indmax,atoms%ntype))
48 :
49 256 : DO i=1,3
50 414 : DO iType=1, atoms%ntype
51 414 : CALL lattHarmsRepToSphHarms(sym, atoms, sphhar, iType, bxc(i)%mt(:,:,iType,1), flm(:,:,iType))
52 : END DO
53 320 : IF (i==1) THEN
54 64 : CALL gradYlm(atoms,flm,grsflm1)
55 256 : ELSE IF (i==2) THEN
56 64 : CALL gradYlm(atoms,flm,grsflm2)
57 : ELSE
58 64 : CALL gradYlm(atoms,flm,grsflm3)
59 : END IF
60 : END DO
61 :
62 64 : DEALLOCATE(flm)
63 :
64 64 : CALL divYlm(grsflm1(:,:indmax,:,:),grsflm2(:,:indmax,:,:),grsflm3(:,:indmax,:,:), divflm)
65 :
66 138 : DO iType=1, atoms%ntype
67 138 : CALL sphHarmsRepToLattHarms(sym, atoms, sphhar, iType, divflm(:,1:indmax,iType), div%mt(:,0:,iType,1))
68 : END DO
69 :
70 64 : DEALLOCATE(divflm,grsflm1,grsflm2,grsflm3)
71 :
72 64 : CALL timestop("MT divergence")
73 :
74 64 : CALL timestart("PW divergence")
75 :
76 312118 : div%pw(:,1)=CMPLX(0.0,0.0)
77 :
78 256 : DO i=1,3
79 936418 : div%pw(:,1)=div%pw(:,1)+ImagUnit*(cell%bmat(i,1)*stars%kv3(1,:)+cell%bmat(i,2)*stars%kv3(2,:)+cell%bmat(i,3)*stars%kv3(3,:))*bxc(i)%pw(:,1)
80 : END DO
81 :
82 64 : CALL timestop("PW divergence")
83 :
84 64 : IF (input%film) THEN
85 0 : CALL timestart("Vac divergence")
86 : !div%vacxy=CMPLX(0.0,0.0)
87 : !div%vacz=0.0
88 0 : CALL vac_grad(vacuum,stars,cell,bxc(1),grad,9*stars%mx1*stars%mx2)
89 : !div%vacxy=div%vacxy+grad(1)%vacxy
90 : !div%vacz=div%vacz+grad(1)%vacz
91 0 : div%vac(:vacuum%nmzxyd,2:,:,:)=div%vac(:vacuum%nmzxyd,2:,:,:)+grad(1)%vac(:vacuum%nmzxyd,2:,:,:)
92 0 : div%vac(:,1,:,:)=div%vac(:,1,:,:)+grad(1)%vac(:,1,:,:)
93 0 : CALL vac_grad(vacuum,stars,cell,bxc(2),grad,9*stars%mx1*stars%mx2)
94 : !div%vacxy=div%vacxy+grad(2)%vacxy
95 : !div%vacz=div%vacz+grad(2)%vacz
96 0 : div%vac(:vacuum%nmzxyd,2:,:,:)=div%vac(:vacuum%nmzxyd,2:,:,:)+grad(2)%vac(:vacuum%nmzxyd,2:,:,:)
97 0 : div%vac(:,1,:,:)=div%vac(:,1,:,:)+grad(2)%vac(:,1,:,:)
98 0 : CALL vac_grad(vacuum,stars,cell,bxc(3),grad,9*stars%mx1*stars%mx2)
99 : !div%vacxy=div%vacxy+grad(3)%vacxy
100 : !div%vacz=div%vacz+grad(3)%vacz
101 0 : div%vac(:vacuum%nmzxyd,2:,:,:)=div%vac(:vacuum%nmzxyd,2:,:,:)+grad(3)%vac(:vacuum%nmzxyd,2:,:,:)
102 0 : div%vac(:,1,:,:)=div%vac(:,1,:,:)+grad(3)%vac(:,1,:,:)
103 0 : CALL timestop("Vac divergence")
104 : END IF
105 :
106 :
107 3072 : END SUBROUTINE divergence
108 :
109 0 : SUBROUTINE vac_grad(vacuum,stars,cell,den,grad,ifftd2)
110 :
111 : USE m_constants
112 : USE m_grdchlh
113 : USE m_fft2d
114 : USE m_types
115 :
116 : IMPLICIT NONE
117 : TYPE(t_vacuum),INTENT(IN) :: vacuum
118 : TYPE(t_stars),INTENT(IN) :: stars
119 : TYPE(t_cell),INTENT(IN) :: cell
120 : TYPE(t_potden),INTENT(IN) :: den
121 : TYPE(t_potden),INTENT(INOUT),DIMENSION(3) :: grad
122 : ! ..
123 : ! .. Scalar Arguments ..
124 : INTEGER, INTENT (IN) :: ifftd2
125 :
126 : ! ..
127 : ! .. Local Scalars ..
128 : INTEGER :: js,nt,i,iq,irec2,nmz0,nmzdiff,ivac,ip
129 : REAL :: rhti,zro,d_15
130 : ! ..
131 : ! .. Local Arrays ..
132 : REAL :: fgz(3)
133 0 : REAL, ALLOCATABLE :: af2(:),bf2(:)
134 0 : REAL, ALLOCATABLE :: rhdx(:),rhdy(:),rhdz(:)
135 0 : REAL, ALLOCATABLE :: rhtdz(:),rhtdzz(:)
136 0 : REAL, ALLOCATABLE :: rxydzr(:),rxydzi(:)
137 0 : REAL, ALLOCATABLE :: rxydzzr(:),rxydzzi(:),rhtxyr(:),rhtxyi(:)
138 0 : REAL, ALLOCATABLE :: rhtxc(:,:),dummy(:)
139 0 : COMPLEX, ALLOCATABLE :: fgxy(:,:),fg(:,:),rxydz(:,:),rxydzz(:),cqpw(:)
140 0 : COMPLEX, ALLOCATABLE :: rdz(:,:), rdzz(:,:)
141 :
142 0 : d_15 = 1.e-15
143 0 : zro = 0.0
144 0 : nt = ifftd2
145 :
146 0 : ALLOCATE (rxydz(vacuum%nmzxy,stars%ng2-1))
147 0 : ALLOCATE (rhtdz(vacuum%nmzd),rhtdzz(vacuum%nmzd))
148 0 : ALLOCATE (rdz(vacuum%nmzd,stars%ng2),rdzz(vacuum%nmzd,stars%ng2))
149 :
150 0 : DO ivac=1,vacuum%nvac
151 :
152 : ! the charge density in vacuum is expanded in 2-dim stars on a mesh
153 : ! in z-direction. the g||.ne.zero-components expand from 1 to nmzxy
154 : ! the g||.eq.zero-components expand from 1 to nmz
155 : ! first we calculate vxc in the warping region
156 :
157 : !
158 : ! calculate first (rhtdz) & second (rhtdzz) derivative of den%vacz(1:nmz)
159 : !
160 :
161 : CALL grdchlh(vacuum%delz,REAL(den%vac(1:vacuum%nmz,1,ivac,1)),&
162 0 : rhtdz(:),rhtdzz)
163 0 : ALLOCATE ( rhtxyr(vacuum%nmzxy), rhtxyi(vacuum%nmzxy),dummy(vacuum%nmzxy) )
164 0 : ALLOCATE ( rxydzr(vacuum%nmzxy), rxydzi(vacuum%nmzxy) )
165 0 : ALLOCATE ( rxydzzr(vacuum%nmzxy),rxydzzi(vacuum%nmzxy) )
166 :
167 0 : DO iq = 2, stars%ng2
168 : !
169 : ! calculate first (rxydz) & second (rxydzz) derivative of den%vacxy:
170 : !
171 0 : DO ip=1,vacuum%nmzxy
172 0 : rhtxyr(ip)=den%vac(ip,iq,ivac,1)
173 : ENDDO
174 0 : CALL grdchlh(vacuum%delz,rhtxyr(:vacuum%nmzxy), rxydzr,rxydzzr)
175 :
176 0 : DO ip=1,vacuum%nmzxy
177 0 : rhtxyi(ip)=aimag(den%vac(ip,iq,ivac,js))
178 : ENDDO
179 :
180 0 : CALL grdchlh(vacuum%delz,rhtxyi(:vacuum%nmzxy), rxydzi,rxydzzi)
181 :
182 0 : DO ip=1,vacuum%nmzxy
183 0 : rdz(ip,iq)=cmplx(rxydzr(ip),rxydzi(ip))
184 : ENDDO
185 :
186 : ENDDO ! loop over 2D stars (iq)
187 :
188 0 : DEALLOCATE ( rhtxyr,rhtxyi,rxydzr,rxydzi,rxydzzr,rxydzzi )
189 0 : DEALLOCATE ( dummy )
190 :
191 0 : ALLOCATE ( rhdx(0:ifftd2-1),rhdy(0:ifftd2-1) )
192 0 : ALLOCATE ( rhdz(0:ifftd2-1))
193 :
194 0 : ALLOCATE ( cqpw(stars%ng2),af2(0:ifftd2-1) )
195 0 : ALLOCATE ( fg(stars%ng2,3),bf2(0:ifftd2-1) )
196 :
197 0 : af2=0.0
198 0 : cqpw = CMPLX(0.0,0.0)
199 0 : DO ip = 1,vacuum%nmzxy
200 : ! loop over warping region
201 :
202 : ! Transform charge and magnetization to real-space.
203 :
204 0 : CALL fft2d(stars, af2(0),bf2, den%vac(ip,:,ivac,1),+1)
205 :
206 : ! calculate derivatives with respect to x,y in g-space
207 : ! and transform them to real-space.
208 :
209 0 : DO iq=2,stars%ng2
210 0 : cqpw(iq)=ImagUnit*den%vac(ip,iq,ivac,js)
211 : ENDDO
212 :
213 : ! d(rho)/atoms%dx is obtained by a FFT of i*gx*den%vac
214 : ! dn/x = FFT(i*gx*den%vac)
215 0 : rhti = 0.0
216 0 : CALL fft2d(stars, rhdx(0),bf2, cqpw,+1,firstderiv=[1.,0.0,0.],cell=cell)
217 :
218 : ! dn/dy = FFT(i*gy*den%vac)
219 0 : rhti = 0.0
220 : CALL fft2d( &
221 0 : stars, rhdy(0),bf2, cqpw, +1,firstderiv=[0.,1.0,0.],cell=cell)
222 :
223 : ! dn/dz = FFT(rdz)
224 0 : rhti = 0.0
225 0 : CALL fft2d(stars, rhdz(0),bf2, rdz(ip,:), +1)
226 :
227 : ! set minimal value of af2 to 1.0e-15
228 : ! af2=max(af2,10e-13)
229 0 : WHERE (af2<d_15) af2=d_15
230 :
231 : ! ----> 2-d back fft to g space
232 0 : bf2=0.0
233 0 : CALL fft2d(stars, rhdx,bf2, fg(:,1), -1)
234 0 : CALL fft2d(stars, rhdy,bf2, fg(:,2), -1)
235 0 : CALL fft2d(stars, rhdz,bf2, fg(:,3), -1)
236 :
237 : ! All the components are added to grad%vac
238 0 : DO irec2 = 1,stars%ng2
239 0 : grad(1)%vac(ip,irec2,ivac,1)=grad(1)%vac(ip,irec2,ivac,1)+fg(irec2,1)
240 0 : grad(2)%vac(ip,irec2,ivac,1)=grad(2)%vac(ip,irec2,ivac,1)+fg(irec2,2)
241 0 : grad(3)%vac(ip,irec2,ivac,1)=grad(3)%vac(ip,irec2,ivac,1)+fg(irec2,3)
242 : ENDDO
243 :
244 : END DO ! ip=1,vacuum%nmzxy
245 0 : DEALLOCATE ( rhdx,rhdy,rhdz)
246 0 : DEALLOCATE ( cqpw,fgxy)
247 :
248 : ! now treat the non-warping region
249 :
250 : nmzdiff = vacuum%nmz - vacuum%nmzxy
251 :
252 : ! The non-warping region runs from nmzxy+1 to nmz.
253 : ! The values from nmz0 to nmzxy are taken into account in order
254 : ! to get the real-space derivative smooth around nmzxy+1.
255 :
256 : nmz0= vacuum%nmzxy+1+(6/2)-6
257 : IF (nmz0 <= 0) THEN ! usually vacuum%nmzxy>6
258 : nmz0= 1
259 : END IF
260 :
261 : DEALLOCATE ( af2)
262 :
263 : DO ip = vacuum%nmzxy + 1,vacuum%nmz
264 : grad(3)%vac(ip,1,ivac,1) = grad(3)%vac(ip,1,ivac,1) + rhtdz(ip-vacuum%nmzxy)
265 : ENDDO
266 :
267 0 : DEALLOCATE ( bf2)
268 :
269 : ENDDO ! loop over vacua (ivac)
270 :
271 0 : END SUBROUTINE vac_grad
272 :
273 2 : SUBROUTINE divpotgrad(input,stars,atoms,sphhar,vacuum,sym,cell,noco,pot,grad)
274 :
275 : USE m_types
276 : USE m_lattHarmsSphHarmsConv
277 : USE m_gradYlm
278 : USE m_constants
279 :
280 : !--------------------------------------------------------------------------
281 : ! Use the interstitial/vacuum gradient subroutine and an external MT-gra-
282 : ! dient routine from juPhon to assemble the gradient of a potenital into a
283 : ! t_potden variable. The MT-gradient is first calculated in sperical har-
284 : ! monics coefficients.
285 : !--------------------------------------------------------------------------
286 :
287 : IMPLICIT NONE
288 :
289 : TYPE(t_input), INTENT(IN) :: input
290 : TYPE(t_stars),INTENT(IN) :: stars
291 : TYPE(t_atoms), INTENT(IN) :: atoms
292 : TYPE(t_sphhar), INTENT(IN) :: sphhar
293 : TYPE(t_vacuum),INTENT(IN) :: vacuum
294 : TYPE(t_sym), INTENT(IN) :: sym
295 : TYPE(t_cell),INTENT(IN) :: cell
296 : TYPE(t_noco), INTENT(IN) :: noco
297 : TYPE(t_potden), INTENT(IN) :: pot
298 : TYPE(t_potden), dimension(3), INTENT(INOUT) :: grad
299 :
300 2 : TYPE(t_potden) :: denloc
301 : INTEGER :: i,iType,indmax,lh,lhmax
302 2 : COMPLEX, ALLOCATABLE :: flm(:,:,:),grsflm(:,:,:,:) ! (iR,lm,n[,x,i])
303 :
304 2 : CALL timestart("MT potential gradient")
305 2 : indmax=(atoms%lmaxd+1)**2
306 :
307 10 : ALLOCATE(flm(atoms%jmtd,indmax,atoms%ntype))
308 :
309 2 : denloc=pot
310 :
311 6 : DO iType=1,atoms%ntype
312 4 : lhmax=sphhar%nlh(sym%ntypsy(atoms%firstAtom(iType)))
313 328 : DO lh=0, lhmax
314 245596 : denloc%mt(:,lh,iType,1) = denloc%mt(:,lh,iType,1)*atoms%rmsh(:, iType)**2
315 : END DO ! lh
316 6 : CALL lattHarmsRepToSphHarms(sym, atoms, sphhar, iType, denloc%mt(:,:,iType,1), flm(:,:,iType))
317 : END DO
318 :
319 : CALL gradYlm(atoms,flm,grsflm)
320 :
321 2 : DEALLOCATE(flm)
322 :
323 8 : DO i=1,3
324 20 : DO iType=1,atoms%ntype
325 736794 : CALL sphHarmsRepToLattHarms(sym, atoms, sphhar, iType, grsflm(:,1:indmax,iType,i)/(4.0*pi_const), grad(i)%mt(:,0:,iType,1))
326 : END DO
327 : END DO
328 :
329 2 : DEALLOCATE(grsflm)
330 :
331 2 : CALL timestop("MT potential gradient")
332 :
333 2 : CALL timestart("PW potential gradient")
334 :
335 8 : DO i=1,3
336 18434 : grad(i)%pw(:,1)=ImagUnit*(cell%bmat(i,1)*stars%kv3(1,:)+cell%bmat(i,2)*stars%kv3(2,:)+cell%bmat(i,3)*stars%kv3(3,:))*pot%pw(:,1)/(4.0*pi_const)
337 : END DO
338 :
339 :
340 2 : CALL timestop("PW potential gradient")
341 :
342 2 : IF (input%film) THEN
343 0 : CALL timestart("Vac potential gradient")
344 0 : CALL vac_grad(vacuum,stars,cell,pot,grad,9*stars%mx1*stars%mx2)
345 0 : CALL timestart("Vac potential gradient")
346 : END IF
347 :
348 4 : END SUBROUTINE divpotgrad
349 : END MODULE m_divergence
|
