Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (C) 2020 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_vac_tofrom_grid
7 : INTEGER,PARAMETER :: fixed_ndvgrd=6
8 :
9 : CONTAINS
10 44 : subroutine vac_to_grid(dograds,ifftd2,jspins,vacuum,l_noco,cell,vacnew,stars,rho,grad,rhoim)
11 :
12 :
13 : !-----------------------------------------------------------------------
14 : ! instead of vvacxcor.f: the different exchange-correlation
15 : ! potentials defined through the key icorr are called through
16 : ! the driver subroutine vxcallg.f, subroutines vectorized
17 : ! in case of total = .true. calculates the ex-corr. energy
18 : ! density through the driver subroutine excallg.f
19 : ! ** r.pentcheva 08.05.96
20 : !-----------------------------------------------------------------------
21 :
22 : USE m_types
23 : use m_constants
24 : USE m_grdrsvac
25 : USE m_grdchlh
26 : USE m_mkgz
27 : USE m_mkgxyz3
28 : !
29 : !
30 : USE m_fft2d
31 :
32 : IMPLICIT NONE
33 : logical,intent(in) :: dograds
34 : INTEGER,INTENT(IN) :: jspins
35 : TYPE(t_vacuum),INTENT(IN) :: vacuum
36 : LOGICAL,INTENT(IN) :: l_noco
37 : TYPE(t_stars),INTENT(IN) :: stars
38 : TYPE(t_cell),INTENT(IN) :: cell
39 : COMPLEX,INTENT(IN) :: vacnew(:,:,:,:)
40 : TYPE(t_gradients),INTENT(INOUT)::grad
41 : real,intent(OUT) :: rho(:,:)
42 : real, optional, allocatable, intent(out) :: rhoim(:,:)
43 : ! .. Scalar Arguments ..
44 : INTEGER, INTENT (IN) :: ifftd2
45 :
46 : ! ..
47 : ! .. Local Scalars ..
48 : INTEGER :: js,nt,i,iq,irec2,nmz0,nmzdiff,ivac,ip,idx,idx1,idx_loc
49 : REAL :: rhti,zro,fgz,rhmnv,d_15,rd
50 : ! ..
51 : ! .. Local Arrays ..
52 44 : REAL, ALLOCATABLE :: bf2(:)
53 44 : REAL, ALLOCATABLE :: rhdx(:,:),rhdy(:,:),rhdz(:,:)
54 44 : REAL, ALLOCATABLE :: rhdxx(:,:),rhdyy(:,:),rhtdz(:,:),rhtdzz(:,:)
55 44 : REAL, ALLOCATABLE :: rhdzz(:,:),rhdyz(:,:),rhdzx(:,:),rhdxy(:,:)
56 44 : REAL, ALLOCATABLE :: rxydzr(:),rxydzi(:)
57 44 : REAL, ALLOCATABLE :: rxydzzr(:),rxydzzi(:),rhtxyr(:),rhtxyi(:)
58 : REAL, ALLOCATABLE :: rhtxc(:,:)
59 44 : COMPLEX, ALLOCATABLE :: rxydz(:,:,:),rxydzz(:,:,:),cqpw(:)
60 44 : COMPLEX, ALLOCATABLE :: rdz(:,:,:),rdzz(:,:,:)
61 :
62 : ! ..
63 : ! for the noco-case only
64 : REAL :: chdens
65 44 : REAL, ALLOCATABLE :: magmom(:,:), dxmagmom(:),ddxmagmom(:,:)
66 44 : REAL, ALLOCATABLE :: dymagmom(:),ddymagmom(:,:), dzmagmom(:,:),ddzmagmom(:,:)
67 44 : REAL, ALLOCATABLE :: mx(:),my(:)
68 : ! .. unused input (needed for other noco GGA-implementations) ..
69 :
70 :
71 :
72 44 : d_15 = 1.e-15
73 44 : zro = 0.0
74 44 : nt = ifftd2
75 44 : idx=1
76 :
77 13870602 : rho = 0.0
78 :
79 132 : ALLOCATE ( bf2(ifftd2) )
80 44 : IF (PRESENT(rhoim)) THEN
81 0 : ALLOCATE(rhoim,mold=rho)
82 0 : rhoim=0.0
83 : ENDIF
84 44 : WRITE (oUnit,'(/'' ifftd2,vacuum%nmz='',2i7)') ifftd2,vacuum%nmz
85 44 : WRITE (oUnit,'('' 9990nmzxy='',2i5)') vacuum%nmzxy
86 :
87 396 : ALLOCATE ( rxydz(vacuum%nmzxy,stars%ng2,jspins),rxydzz(vacuum%nmzxyd,stars%ng2,jspins) )
88 264 : ALLOCATE ( rhtdz(vacuum%nmzd,jspins),rhtdzz(vacuum%nmzd,jspins) )
89 352 : ALLOCATE ( rdz(vacuum%nmzd,stars%ng2,jspins),rdzz(vacuum%nmzd,stars%ng2,jspins))
90 : !ALLOCATE ( fgxy(stars%ng2-1) )
91 2003033 : rxydz = CMPLX(0.0,0.0)
92 2003033 : rxydzz= CMPLX(0.0,0.0)
93 :
94 44 : IF (l_noco) THEN
95 0 : ALLOCATE ( magmom(0:ifftd2-1,vacuum%nmzxy) )
96 0 : ALLOCATE ( dzmagmom(0:ifftd2-1,vacuum%nmzxy) )
97 0 : ALLOCATE ( ddzmagmom(0:ifftd2-1,vacuum%nmzxy) )
98 0 : ALLOCATE ( mx(0:ifftd2-1),my(0:ifftd2-1) )
99 : ENDIF
100 44 : IF ( l_noco .OR. dograds ) THEN
101 102 : ALLOCATE ( rhtxyr(vacuum%nmzxy) )
102 102 : ALLOCATE ( rxydzr(vacuum%nmzxy),rxydzzr(vacuum%nmzxy) )
103 : ENDIF
104 44 : IF (dograds) THEN
105 102 : ALLOCATE ( rhtxyi(vacuum%nmzxy) )
106 68 : ALLOCATE ( rxydzi(vacuum%nmzxy) )
107 68 : ALLOCATE ( rxydzzi(vacuum%nmzxy) )
108 : ENDIF
109 98 : DO ivac=1,vacuum%nvac
110 :
111 : ! the charge density in vacuum is expanded in 2-dim stars on a mesh
112 : ! in z-direction. the g||.ne.zero-components expand from 1 to nmzxy
113 : ! the g||.eq.zero-components expand from 1 to nmz
114 : ! first we calculate vxc in the warping region
115 :
116 :
117 : ! Transform charge and magnetization to real-space.
118 : ! In the collinear case that is done later within
119 : ! another loop over the vacuum-layers in order to
120 : ! save memory.
121 :
122 : idx1=idx
123 : !idx1=(ivac-1)* ( vacuum%nmzxy * ifftd2 + nmzdiff ) + 1
124 5454 : DO ip=1,vacuum%nmzxy
125 13200 : DO js=1,jspins
126 13200 : IF (.NOT.PRESENT(rhoim)) THEN
127 7800 : CALL fft2d(stars, rho(idx1:idx1+ifftd2-1,js),bf2, vacnew(ip,:,ivac,js),+1)
128 : ELSE
129 0 : CALL fft2d(stars, rho(idx1:idx1+ifftd2-1,js),rhoim(idx1:idx1+ifftd2-1,js), vacnew(ip,:,ivac,js),+1)
130 : END IF
131 : END DO
132 5400 : IF (l_noco) THEN
133 0 : CALL fft2d(stars, mx,my, vacnew(ip,:,ivac,3),+1)
134 :
135 0 : DO i=0,ifftd2-1
136 0 : magmom(i,ip)= mx(i)**2 + my(i)**2 + ((rho(i+idx1,1)-rho(i+idx1,2))/2.)**2
137 0 : magmom(i,ip)= SQRT(magmom(i,ip))
138 0 : chdens= rho(i+idx1,1)/2.+rho(i+idx1,2)/2.
139 0 : rho(i+idx1,1)= chdens + magmom(i,ip)
140 0 : rho(i+idx1,2)= chdens - magmom(i,ip)
141 : END DO
142 : ENDIF
143 5454 : idx1=idx1+ifftd2
144 : END DO ! ip=1,vacuum%nmzxy
145 :
146 : ! ENDDO ! ivac
147 : ! DO ivac = 1,nvac
148 :
149 54 : IF (dograds) THEN
150 72 : DO js=1,jspins
151 : !
152 : ! calculate first (rhtdz) & second (rhtdzz) derivative of vacz(1:nmz)
153 : CALL grdchlh(vacuum%delz,REAL(vacnew(1:vacuum%nmz,1,ivac,js)),&
154 9576 : rhtdz(1:,js),rhtdzz(1:,js))
155 9538 : rdz(:,1,js) = rhtdz(1:,js)
156 9538 : rdzz(:,1,js) = rhtdzz(1:,js)
157 16511 : DO iq = 1, stars%ng2-1
158 : !
159 : ! calculate first (rxydz) & second (rxydzz) derivative of vacxy:
160 : !
161 1663773 : DO ip=1,vacuum%nmzxy
162 1663773 : rhtxyr(ip)=real(vacnew(ip,iq+1,ivac,js))
163 : ENDDO
164 16473 : CALL grdchlh(vacuum%delz,rhtxyr(:vacuum%nmzxy), rxydzr,rxydzzr)
165 :
166 1663773 : DO ip=1,vacuum%nmzxy
167 1663773 : rhtxyi(ip)=aimag(vacnew(ip,iq+1,ivac,js))
168 : ENDDO
169 16473 : CALL grdchlh(vacuum%delz,rhtxyi(:vacuum%nmzxy), rxydzi,rxydzzi)
170 :
171 1663811 : DO ip=1,vacuum%nmzxy
172 1647300 : rxydz(ip,iq+1,js)=cmplx(rxydzr(ip),rxydzi(ip))
173 1663773 : rxydzz(ip,iq+1,js)=cmplx(rxydzzr(ip),rxydzzi(ip))
174 : ENDDO
175 :
176 : ENDDO ! loop over 2D stars (iq)
177 1663811 : rdz(:vacuum%nmzxy,2:,js) = rxydz(:vacuum%nmzxy,2:,js)
178 1663845 : rdzz(:vacuum%nmzxy,2:,js) = rxydzz(:vacuum%nmzxy,2:,js)
179 :
180 : ENDDO ! jspins
181 :
182 34 : IF (l_noco) THEN
183 : ! calculate dzmagmom = d magmom / d z and ddzmagmom= d dmagmom / d z
184 :
185 0 : DO i=0,ifftd2-1
186 0 : DO ip=1,vacuum%nmzxy
187 0 : rhtxyr(ip)=magmom(i,ip)
188 : ENDDO
189 0 : CALL grdchlh(vacuum%delz,rhtxyr(1:vacuum%nmzxy), rxydzr,rxydzzr)
190 0 : DO ip=1,vacuum%nmzxy
191 0 : dzmagmom(i,ip)= rxydzr(ip)
192 0 : ddzmagmom(i,ip)= rxydzzr(ip)
193 : ENDDO
194 : END DO
195 : END IF ! l_noco
196 :
197 : ENDIF ! xcpot%igrd.GT.0
198 :
199 : ! WRITE(oUnit,'('' 9990nmzxy='',2i5)') nmzxy
200 :
201 54 : CALL timestart("warp")
202 54 : rd = 0.0
203 : !$OMP PARALLEL DEFAULT(none) &
204 : !$OMP SHARED(vacuum,dograds,jspins,stars,ivac,zro,cell,magmom,vacnew) &
205 : !$OMP SHARED(rhtdz,rhtdzz,rdz,rdzz,rxydz,rxydzz,l_noco,dzmagmom,ddzmagmom,idx) &
206 : !$OMP SHARED(ifftd2,rho,rhoim,grad) &
207 : !$OMP PRIVATE(ip,js,iq,cqpw,bf2,rhti,rhdx,rhdy,rhdz,rhdxx,rhdyy,rhdzz) &
208 : !$OMP PRIVATE(rhdxy,rhdzx,rhdyz,dxmagmom,dymagmom,ddxmagmom,ddymagmom) &
209 54 : !$OMP PRIVATE(chdens,idx_loc)
210 : ALLOCATE ( rhdx(0:ifftd2-1,jspins),rhdy(0:ifftd2-1,jspins) )
211 : ALLOCATE ( rhdz(0:ifftd2-1,jspins),rhdxx(0:ifftd2-1,jspins) )
212 : ALLOCATE ( rhdyy(0:ifftd2-1,jspins),rhdzz(0:ifftd2-1,jspins) )
213 : ALLOCATE ( rhdyz(0:ifftd2-1,jspins),rhdzx(0:ifftd2-1,jspins) )
214 : ALLOCATE ( rhdxy(0:ifftd2-1,jspins))
215 : ALLOCATE ( cqpw(stars%ng2))
216 : IF (l_noco) THEN
217 : ALLOCATE ( dxmagmom(0:ifftd2-1),dymagmom(0:ifftd2-1) )
218 : ALLOCATE ( ddxmagmom(0:ifftd2-1,2),ddymagmom(0:ifftd2-1,2) )
219 : ENDIF
220 : !$OMP DO
221 : DO ip = 1,vacuum%nmzxy
222 : ! loop over warping region
223 :
224 :
225 : IF (dograds) THEN
226 : ! calculate derivatives with respect to x,y in g-space
227 : ! and transform them to real-space.
228 :
229 : DO js = 1,jspins
230 :
231 : DO iq=2,stars%ng2
232 : cqpw(iq)=ImagUnit*vacnew(ip,iq,ivac,js)
233 : ENDDO
234 :
235 : cqpw(1) = CMPLX(0.0,0.0)
236 : ! d(rho)/atoms%dx is obtained by a FFT of i*gx*vacxy
237 : ! (vacz is set to zero and gx is included in
238 : ! dn/atoms = FFT(0,i*gx*vacxy)
239 :
240 :
241 :
242 : CALL fft2d(stars, rhdx(0,js),bf2, cqpw, +1,firstderiv=[1.,0.,0.],cell=cell)
243 : !TODO & pgft2x)
244 :
245 : cqpw(1) = CMPLX(0.0,0.0)
246 : CALL fft2d(stars, rhdy(0,js),bf2, cqpw, +1,firstderiv=[0.,1.,0.],cell=cell) ! dn/dy = FFT(0,i*gy*vacxy)&
247 :
248 : CALL fft2d(stars, rhdz(0,js),bf2, rdz(ip,:,js), +1) ! dn/dz = FFT(rhtdz,rxydz)&
249 :
250 : DO iq=2,stars%ng2
251 : cqpw(iq)=-vacnew(ip,iq,ivac,js)
252 : ENDDO
253 :
254 : cqpw(1) = CMPLX(0.0,0.0)
255 : CALL fft2d(stars, rhdxx(0,js),bf2, cqpw, +1,firstderiv=[1.0,0.,0.],secondderiv=[1.0,0.,0.],cell=cell) ! d2n/dx2 = FFT(0,-gx^2*vacxy)&
256 :
257 : cqpw(1) = CMPLX(0.0,0.0)
258 : CALL fft2d(stars, rhdyy(0,js),bf2, cqpw, +1,firstderiv=[0.,1.0,0.],secondderiv=[0.,1.0,0.],cell=cell) ! d2n/dy2 = FFT(0,-gy^2*vacxy)&
259 :
260 : rhti = 0.0
261 : CALL fft2d(stars, rhdzz(0,js),bf2, rdzz(ip,:,js), +1) ! d2n/dz2 = FFT(rhtdzz,rxydzz)&
262 :
263 :
264 : DO iq=2,stars%ng2
265 : cqpw(iq)=ImagUnit*rxydz(ip,iq,js)
266 : ENDDO
267 :
268 : cqpw(1) = CMPLX(0.0,0.0)
269 : CALL fft2d(stars, rhdyz(0,js),bf2, cqpw, +1,firstderiv=[0.,1.0,0.],cell=cell) ! d2n/dyz = FFT(0,i*gy*rxydz)&
270 :
271 : cqpw(1) = CMPLX(0.0,0.0)
272 : CALL fft2d(stars, rhdzx(0,js),bf2, cqpw, +1,firstderiv=[1.,0.0,0.],cell=cell) ! d2n/dzx = FFT(0,i*gx*rxydz)&
273 :
274 : DO iq=2,stars%ng2
275 : cqpw(iq)=-vacnew(ip,iq,ivac,js)
276 : ENDDO
277 :
278 : cqpw(1) = CMPLX(0.0,0.0)
279 : CALL fft2d(stars, rhdxy(0,js),bf2, cqpw, +1,firstderiv=[0.,1.0,0.],secondderiv=[1.,0.0,0.],cell=cell) ! d2n/dxy = FFT(0,-gx*gy*vacxy)&
280 :
281 : END DO ! js=1,jspins
282 :
283 :
284 : IF (l_noco) THEN
285 : ! ! In non-collinear calculations the derivatives of |m| are calculated
286 : ! ! in real-space. The derivatives of the charge density, that are
287 : ! ! already calculated in g-space, will be used.
288 :
289 : CALL grdrsvac(magmom(0,ip),cell%bmat,3*stars%mx1,3*stars%mx2,fixed_ndvgrd, dxmagmom,dymagmom)
290 : DO i=0,ifftd2-1
291 : chdens= rhdx(i,1)/2.+rhdx(i,2)/2.
292 : rhdx(i,1)= chdens + dxmagmom(i)
293 : rhdx(i,2)= chdens - dxmagmom(i)
294 : chdens= rhdy(i,1)/2.+rhdy(i,2)/2.
295 : rhdy(i,1)= chdens + dymagmom(i)
296 : rhdy(i,2)= chdens - dymagmom(i)
297 : chdens= rhdz(i,1)/2.+rhdz(i,2)/2.
298 : rhdz(i,1)= chdens + dzmagmom(i,ip)
299 : rhdz(i,2)= chdens - dzmagmom(i,ip)
300 : END DO
301 :
302 : CALL grdrsvac(dxmagmom,cell%bmat,3*stars%mx1,3*stars%mx2,fixed_ndvgrd, &
303 : ddxmagmom(0,1),ddymagmom(0,1))
304 : CALL grdrsvac(&
305 : dymagmom,cell%bmat,3*stars%mx1,3*stars%mx2,fixed_ndvgrd,ddxmagmom(0,2),ddymagmom(0,2))
306 : DO i=0,ifftd2-1
307 : chdens= rhdxx(i,1)/2.+rhdxx(i,2)/2.
308 : rhdxx(i,1)= chdens + ddxmagmom(i,1)
309 : rhdxx(i,2)= chdens - ddxmagmom(i,1)
310 : chdens= rhdyy(i,1)/2.+rhdyy(i,2)/2.
311 : rhdyy(i,1)= chdens + ddymagmom(i,2)
312 : rhdyy(i,2)= chdens - ddymagmom(i,2)
313 : chdens= rhdxy(i,1)/2.+rhdxy(i,2)/2.
314 : rhdxy(i,1)= chdens + ( ddxmagmom(i,2) + ddymagmom(i,1) )/2.
315 : rhdxy(i,2)= chdens - ( ddxmagmom(i,2) + ddymagmom(i,1) )/2.
316 : END DO
317 : CALL grdrsvac(dzmagmom(0,ip),cell%bmat,3*stars%mx1,3*stars%mx2,fixed_ndvgrd, &
318 : ddxmagmom(0,1),ddymagmom(0,1))
319 : DO i=0,ifftd2-1
320 : chdens= rhdzx(i,1)/2.+rhdzx(i,2)/2.
321 : rhdzx(i,1)= chdens + ddxmagmom(i,1)
322 : rhdzx(i,2)= chdens - ddxmagmom(i,1)
323 : chdens= rhdyz(i,1)/2.+rhdyz(i,2)/2.
324 : rhdyz(i,1)= chdens + ddymagmom(i,1)
325 : rhdyz(i,2)= chdens - ddymagmom(i,1)
326 : chdens= rhdzz(i,1)/2.+rhdzz(i,2)/2.
327 : rhdzz(i,1)= chdens + ddzmagmom(i,ip)
328 : rhdzz(i,2)= chdens - ddzmagmom(i,ip)
329 : END DO
330 :
331 : END IF ! l_noco
332 : !
333 : idx_loc = idx + (ip-1)* ifftd2
334 : CALL mkgxyz3(rho(idx_loc:idx_loc+ifftd2-1,:),rhdx,rhdy, rhdz,rhdxx,rhdyy,rhdzz,rhdyz,rhdzx,rhdxy, idx_loc-1,grad)
335 : !
336 : END IF ! vxc_is_gga
337 : !
338 : ! set minimal value of af2 to 1.0e-13
339 : !
340 :
341 : ! rho=max(rho,1e-13)
342 :
343 : END DO ! ip=1,vacuum%nmzxy
344 : !$OMP END DO
345 : DEALLOCATE ( rhdx,rhdy )
346 : DEALLOCATE ( rhdz,rhdxx )
347 : DEALLOCATE ( rhdyy,rhdzz )
348 : DEALLOCATE ( rhdyz,rhdzx )
349 : DEALLOCATE ( rhdxy,cqpw )
350 : IF (l_noco) THEN
351 : DEALLOCATE ( dxmagmom,dymagmom )
352 : DEALLOCATE ( ddxmagmom,ddymagmom )
353 : ENDIF
354 : !$OMP END PARALLEL
355 54 : idx = idx + vacuum%nmzxy * ifftd2
356 54 : CALL timestop("warp")
357 :
358 : ! now treat the non-warping region
359 :
360 :
361 : ! The non-warping region runs from nmzxy+1 to nmz.
362 : ! The values from nmz0 to nmzxy are taken into account in order
363 : ! to get the real-space derivative smooth around nmzxy+1.
364 54 : nmz0= max(1,vacuum%nmzxy+1+(fixed_ndvgrd/2)-fixed_ndvgrd)
365 54 : nmzdiff = vacuum%nmz - nmz0+1
366 : ! WRITE(oUnit,'(/'' 9992excz''/(8f15.7))') (excz(ip,1),ip=1,nmz)
367 54 : WRITE(oUnit,'(/'' 9992nmzdiff='',i5)') nmzdiff
368 :
369 :
370 : !idx = (ivac-1)* ( vacuum%nmzxy * ifftd2 + nmzdiff ) + ip*ifftd2 + 1
371 8316 : DO ip=nmz0,vacuum%nmz
372 8316 : IF (.not. l_noco) THEN
373 20196 : DO js=1,jspins
374 : !rho(idx+ip-nmz0,js)= REAL(vacz(ip,ivac,js))
375 11934 : rho(idx+ip-nmz0,js)= REAL(vacnew(ip,1,ivac,js))
376 20196 : IF (PRESENT(rhoim)) rhoim(idx+ip-nmz0,js)= AIMAG(vacnew(ip,1,ivac,js))
377 : END DO
378 : ELSE
379 : !mx(0) = REAL(vacz(ip,ivac,3))
380 0 : mx(0) = REAL(vacnew(ip,1,ivac,3))
381 : !my(0) = AIMAG(vacz(ip,ivac,3))
382 0 : my(0) = AIMAG(vacnew(ip,1,ivac,3))
383 : !chdens= (REAL(vacz(ip,ivac,1))+REAL(vacz(ip,ivac,2)))/2
384 0 : chdens= (REAL(vacnew(ip,1,ivac,1))+REAL(vacnew(ip,1,ivac,2)))/2
385 : !magmom(0,1)= mx(0)**2 + my(0)**2 + ((REAL(vacz(ip,ivac,1))-REAL(vacz(ip,ivac,2)))/2.)**2
386 0 : magmom(0,1)= mx(0)**2 + my(0)**2 + ((REAL(vacnew(ip,1,ivac,1))-REAL(vacnew(ip,1,ivac,2)))/2.)**2
387 0 : magmom(0,1)= SQRT(magmom(0,1))
388 0 : rho(idx+ip-nmz0,1)= chdens + magmom(0,1)
389 0 : rho(idx+ip-nmz0,2)= chdens - magmom(0,1)
390 : END IF
391 : END DO
392 54 : IF (dograds) THEN
393 34 : IF (l_noco) THEN
394 0 : DO js=1,jspins
395 0 : CALL grdchlh(vacuum%delz,rho(idx:idx+nmzdiff-1,js),rhtdz(nmz0:,js),rhtdzz(nmz0:,js))
396 : END DO
397 :
398 : END IF
399 :
400 : ! calculate the quantities such as abs(grad(rho)),.. used in
401 : !c evaluating the gradient contributions to potential and
402 : !c energy.
403 :
404 :
405 :
406 : CALL mkgz(nmzdiff,jspins, rho(nmz0:,1),rho(nmz0:,jspins),&
407 : rhtdz(nmz0:,1),rhtdz(nmz0:,jspins),rhtdzz(nmz0:,1),&
408 34 : rhtdzz(nmz0:,jspins),idx-1,grad)
409 :
410 : ENDIF
411 :
412 : ! calculate vxc for z now beyond warping region
413 98 : idx=idx+nmzdiff
414 : ENDDO ! loop over vacua (ivac)
415 :
416 :
417 :
418 44 : END SUBROUTINE vac_to_grid
419 :
420 88 : subroutine vac_from_grid(stars,vacuum,v_xc,ifft2d,vac)
421 : use m_types_stars
422 : use m_types_vacuum
423 : use m_fft2d
424 : type(t_stars),intent(in) :: stars
425 : type(t_vacuum),intent(in) :: vacuum
426 :
427 : real, INTENT(IN) :: v_xc(:,:)
428 : INTEGER,INTENT(IN) :: ifft2d
429 : complex,intent(INOUT) :: vac(:,:,:,:)
430 :
431 88 : COMPLEX, ALLOCATABLE :: fg(:)
432 88 : REAL, ALLOCATABLE :: bf2(:)
433 : INTEGER :: js,irec2,idx,ivac,ip,nmz0
434 :
435 440 : ALLOCATE ( fg(stars%ng2),bf2(ifft2d) )
436 :
437 190 : DO js = 1,size(v_xc,2)
438 102 : idx=1
439 322 : DO ivac=1,vacuum%nvac
440 13332 : DO ip=1,vacuum%nmzxy
441 : !
442 : ! ----> 2-d back fft to g space
443 : !
444 25249200 : bf2=0.0
445 13200 : CALL fft2d(stars, v_xc(idx:idx-1+ifft2d,js),bf2, fg, -1)
446 13200 : idx=idx+ifft2d
447 : ! ----> and add vxc to coulomb potential
448 : ! the g||.eq.zero component is added to vxc%vacz
449 : !
450 13200 : vac(ip,1,ivac,js) = fg(1) + vac(ip,1,ivac,js)
451 : !
452 : ! the g||.ne.zero components are added to vxc%vacxy
453 : !
454 4218932 : DO irec2 = 2,stars%ng2
455 4218800 : vac(ip,irec2,ivac,js)=vac(ip,irec2,ivac,js)+fg(irec2)
456 : ENDDO
457 : enddo
458 :
459 132 : nmz0= max(1,vacuum%nmzxy+1+(fixed_ndvgrd/2)-fixed_ndvgrd)
460 :
461 20430 : DO ip = nmz0,vacuum%nmz
462 20196 : if (ip>vacuum%nmzxy) vac(ip,1,ivac,js) = vac(ip,1,ivac,js) + v_xc(idx,js)
463 20328 : idx=idx+1
464 : ENDDO
465 : END DO
466 : ENDDO
467 :
468 88 : end subroutine
469 : end module
|
