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 :
7 : MODULE m_pwden
8 : CONTAINS
9 7322 : SUBROUTINE pwden(stars, kpts, banddos, input, fmpi, noco, nococonv, cell, atoms, sym, &
10 7322 : ikpt, jspin, lapw, ne, ev_list, we, eig, den, results, f_b8, zMat, dos, q_dfpt, lapwq, we1, zMat1, iDir, &
11 : lapwmq,zMat1m)
12 : !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
13 : ! In this subroutine the star function expansion coefficients of
14 : ! the plane wave charge density is determined.
15 : !
16 : ! This subroutine is called for each k-point and each spin.
17 : !
18 : !
19 : ! INPUT: eigen vectors
20 : ! reciprocal lattice information
21 : ! Brillouine zone sampling
22 : ! FFT information
23 : !
24 : ! OUTPUT: den%pw(s)
25 : !
26 : ! Stefan Bl"ugel, JRCAT, Feb. 1997
27 : ! Gustav Bihlmayer, UniWien
28 : !
29 : ! In non-collinear calculations the density becomes a hermitian 2x2
30 : ! matrix. This subroutine generates this density matrix in the
31 : ! interstitial region. The diagonal elements of this matrix
32 : ! (n_11 & n_22) are stored in den%pw, while the real and imaginary part
33 : ! of the off-diagonal element are store in den%pw(:,3).
34 : !
35 : ! Philipp Kurz 99/07
36 : !
37 : ! Subtroutine was refactored in 2020. GM.
38 : !
39 : !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
40 :
41 : !DEC$ NOOPTIMIZE
42 : USE m_types
43 : USE m_types_dos
44 : USE m_constants
45 : USE m_forceb8
46 : USE m_pwint
47 : USE m_juDFT
48 : USE m_types_fftGrid
49 : USE m_fft_interface
50 :
51 : IMPLICIT NONE
52 :
53 : TYPE(t_lapw), INTENT(IN) :: lapw
54 : TYPE(t_mpi), INTENT(IN) :: fmpi
55 :
56 : TYPE(t_banddos), INTENT(IN) :: banddos
57 : TYPE(t_input), INTENT(IN) :: input
58 : TYPE(t_noco), INTENT(IN) :: noco
59 : TYPE(t_nococonv),INTENT(IN) :: nococonv
60 : TYPE(t_sym), INTENT(IN) :: sym
61 : TYPE(t_stars), INTENT(IN) :: stars
62 : TYPE(t_cell), INTENT(IN) :: cell
63 : TYPE(t_kpts), INTENT(IN) :: kpts
64 : TYPE(t_atoms), INTENT(IN) :: atoms
65 : TYPE(t_mat), INTENT(IN) :: zMat
66 : TYPE(t_potden), INTENT(INOUT) :: den
67 : TYPE(t_results), INTENT(INOUT) :: results
68 : TYPE(t_dos), INTENT(INOUT) :: dos
69 :
70 : REAL, INTENT(IN) :: we(:) !(nobd)
71 : REAL, INTENT(IN) :: eig(:) !(input%neig)
72 : INTEGER, INTENT(IN) :: ne
73 : INTEGER, INTENT(IN) :: ev_list(ne)
74 : INTEGER, INTENT(IN) :: ikpt
75 : INTEGER, INTENT(IN) :: jspin
76 : COMPLEX, INTENT(INOUT) :: f_b8(3, atoms%ntype)
77 :
78 : REAL, OPTIONAL, INTENT(IN) :: q_dfpt(3), we1(:)
79 : TYPE(t_mat), OPTIONAL, INTENT(IN) :: zMat1
80 : TYPE(t_lapw), OPTIONAL, INTENT(IN) :: lapwq
81 : INTEGER, OPTIONAL, INTENT(IN) :: iDir
82 :
83 : TYPE(t_lapw), OPTIONAL, INTENT(IN) :: lapwmq
84 : TYPE(t_mat), OPTIONAL, INTENT(IN) :: zMat1m
85 :
86 : ! local variables
87 527184 : TYPE(t_fftGrid) :: state, stateB, stateq, stateBq, StateDeriv, ekinGrid, chargeDen, rhomatGrid(4), statemq, stateBmq
88 51254 : TYPE(t_fftGrid) :: stepFct
89 : INTEGER nu, iv, ir, istr, i, j
90 : INTEGER idens, ndens, ispin, iGp, iG, gVec(3), gInd
91 : REAL q0, q0_11, q0_22, norm, xk(3)
92 : REAL s, stateRadius, stateFFTRadius, stateFFTExtendedRadius
93 : COMPLEX x
94 : REAL, PARAMETER:: tol_3 = 1.0e-3
95 : LOGICAL forw, l_dfpt, l_minusq
96 :
97 : ! local arrays
98 7322 : INTEGER, ALLOCATABLE :: stateIndices(:)
99 7322 : INTEGER, ALLOCATABLE :: stateBIndices(:)
100 7322 : INTEGER, ALLOCATABLE :: stateqIndices(:)
101 7322 : INTEGER, ALLOCATABLE :: stateBqIndices(:)
102 7322 : INTEGER, ALLOCATABLE :: statemqIndices(:)
103 : INTEGER, ALLOCATABLE :: stateBmqIndices(:)
104 7322 : INTEGER, ALLOCATABLE :: fieldSphereIndices(:)
105 7322 : REAL wtf(ne), wtf1(ne)
106 7322 : COMPLEX tempState(lapw%nv(jspin)), starCharges(stars%ng3), z0(lapw%nv(jspin))
107 7322 : COMPLEX, ALLOCATABLE :: cwk(:), ecwk(:)
108 :
109 : ! subroutines
110 : REAL dznrm2, dnrm2
111 :
112 : !-------> ABBREVIATIONS
113 : !
114 : ! ne : number of occupied states
115 : ! nv : number of g-components in eigenstate
116 : ! cwk : complex work array: charge density in g-space (as stars)
117 : ! den%pw : charge density stored as stars
118 : ! we : weights for the BZ-integration for a particular k-point
119 : ! omtil : volume (slab) unit cell, between -.5*D_tilde and +.5*D_tilde
120 : ! nstr : number of members (arms) of reciprocal lattice (g) vector
121 : ! of each star.
122 : ! ng3_fft: number of stars in the charge density FFT-box
123 : ! ng3 : number of 3 dim. stars in the charge density sphere defined
124 : ! by gmax
125 :
126 7322 : CALL timestart("pwden")
127 :
128 7322 : l_dfpt = PRESENT(q_dfpt)
129 :
130 7322 : l_minusq = PRESENT(lapwmq)
131 :
132 1622775 : stateRadius = MAXVAL(ABS(lapw%rk(:,:)))
133 29288 : stateRadius = stateRadius + SQRT(DOT_PRODUCT(kpts%bk(:,ikpt),kpts%bk(:,ikpt)))
134 9410 : IF (noco%l_noco) stateRadius = stateRadius + SQRT(DOT_PRODUCT(nococonv%qss(:),nococonv%qss(:)))
135 7322 : IF (l_dfpt) stateRadius = stateRadius + SQRT(DOT_PRODUCT(q_dfpt,q_dfpt))
136 :
137 7322 : stateFFTRadius = 2.0*stateRadius
138 7322 : stateFFTExtendedRadius = 2.0*stateRadius
139 :
140 7322 : IF (noco%l_noco.OR.noco%l_soc) THEN
141 1758 : IF (banddos%dos .OR. banddos%vacdos .OR. input%cdinf.OR.banddos%band) THEN
142 2 : stateFFTExtendedRadius = 3.0*stateRadius+0.1
143 2 : CALL stepFct%init(cell,sym,stateFFTExtendedRadius+0.001)
144 2 : CALL stepFct%putFieldOnGrid(stars, stars%ustep, cell,stateFFTRadius+0.0005)
145 2 : CALL stepFct%fillFieldSphereIndexArray(stars, stateFFTRadius+0.0008, fieldSphereIndices)
146 2 : CALL fft_interface(3, stepFct%dimensions(:), stepFct%grid, .FALSE., fieldSphereIndices)
147 : END IF
148 : END IF
149 :
150 29288 : ALLOCATE (cwk(stars%ng3), ecwk(stars%ng3))
151 :
152 7322 : IF (noco%l_noco) THEN
153 696 : CALL state%init(cell,sym,stateFFTExtendedRadius+0.001)
154 696 : CALL stateB%init(cell,sym,stateFFTExtendedRadius+0.001)
155 696 : IF (l_dfpt) THEN
156 0 : CALL stateq%init(cell,sym,stateFFTExtendedRadius+0.001)
157 0 : CALL stateBq%init(cell,sym,stateFFTExtendedRadius+0.001)
158 0 : IF (l_minusq) THEN
159 0 : CALL statemq%init(cell,sym,stateFFTExtendedRadius+0.001)
160 0 : CALL stateBmq%init(cell,sym,stateFFTExtendedRadius+0.001)
161 : END IF
162 : END IF
163 3480 : DO i = 1, 4
164 2784 : CALL rhomatGrid(i)%init(cell,sym,stateFFTExtendedRadius+0.001)
165 13288776 : rhomatGrid(i)%grid(:) = CMPLX(0.0,0.0)
166 : END DO
167 696 : CALL rhomatGrid(1)%fillFieldSphereIndexArray(stars, stateFFTRadius+0.0008, fieldSphereIndices)
168 696 : IF (noco%l_ss) THEN
169 102 : ALLOCATE(stateIndices(lapw%nv(1)))
170 102 : ALLOCATE(stateBIndices(lapw%nv(2)))
171 34 : CALL state%fillStateIndexArray(lapw,1,stateIndices)
172 34 : CALL stateB%fillStateIndexArray(lapw,2,stateBIndices)
173 34 : IF (l_dfpt) THEN
174 0 : ALLOCATE(stateqIndices(lapwq%nv(1)))
175 0 : ALLOCATE(stateBqIndices(lapwq%nv(2)))
176 0 : CALL stateq%fillStateIndexArray(lapwq,1,stateqIndices)
177 0 : CALL stateBq%fillStateIndexArray(lapwq,2,stateBqIndices)
178 : END IF
179 : ELSE
180 1986 : ALLOCATE(stateIndices(lapw%nv(jspin)))
181 1324 : ALLOCATE(stateBIndices(lapw%nv(jspin)))
182 662 : CALL state%fillStateIndexArray(lapw,jspin,stateIndices)
183 662 : CALL stateB%fillStateIndexArray(lapw,jspin,stateBIndices)
184 662 : IF (l_dfpt) THEN
185 0 : ALLOCATE(stateqIndices(lapwq%nv(jspin)))
186 0 : ALLOCATE(stateBqIndices(lapwq%nv(jspin)))
187 0 : CALL stateq%fillStateIndexArray(lapwq,jspin,stateqIndices)
188 0 : CALL stateBq%fillStateIndexArray(lapwq,jspin,stateBqIndices)
189 : END IF
190 : ENDIF
191 : ELSE
192 6626 : CALL state%init(cell,sym,stateFFTExtendedRadius+0.001)
193 6626 : IF (l_dfpt) THEN
194 0 : CALL stateq%init(cell,sym,stateFFTExtendedRadius+0.001)
195 0 : IF (l_minusq) THEN
196 0 : CALL statemq%init(cell,sym,stateFFTExtendedRadius+0.001)
197 : END IF
198 : END IF
199 6626 : CALL chargeDen%init(cell,sym,stateFFTExtendedRadius+0.001)
200 22130592 : chargeDen%grid(:) = CMPLX(0.0,0.0)
201 : ! TODO: Shouldn't there be a starsq here for DFPT?
202 6626 : CALL chargeDen%fillFieldSphereIndexArray(stars, stateFFTRadius+0.0008, fieldSphereIndices)
203 6626 : IF (input%l_f) THEN
204 62 : CALL stateDeriv%init(cell,sym,stateFFTExtendedRadius+0.001)
205 62 : CALL ekinGrid%init(cell,sym,stateFFTExtendedRadius+0.001)
206 507187 : ekinGrid%grid(:) = CMPLX(0.0,0.0)
207 : END IF
208 19878 : ALLOCATE(stateIndices(lapw%nv(jspin)))
209 6626 : CALL state%fillStateIndexArray(lapw,jspin,stateIndices)
210 6626 : IF (l_dfpt) THEN
211 0 : ALLOCATE(stateqIndices(lapwq%nv(jspin)))
212 0 : CALL stateq%fillStateIndexArray(lapwq,jspin,stateqIndices)
213 0 : IF (l_minusq) THEN
214 0 : ALLOCATE(statemqIndices(lapwmq%nv(jspin)))
215 0 : CALL statemq%fillStateIndexArray(lapwmq,jspin,statemqIndices)
216 : END IF
217 : END IF
218 : ENDIF
219 :
220 0 : IF (.NOT.l_dfpt) THEN
221 : ! g=0 star: calculate the charge for this k-point and spin
222 : ! analytically to test the quality of FFT
223 7322 : q0 = 0.0
224 7322 : q0_11 = 0.0
225 7322 : q0_22 = 0.0
226 7322 : IF (noco%l_noco) THEN
227 696 : IF (.NOT. zmat%l_real) THEN
228 9505 : DO nu = 1, ne
229 8809 : norm = dznrm2(lapw%nv(1),zMat%data_c(1:, nu),1) ! dznrm2 returns the 2-norm of the vector.
230 8809 : q0_11 = q0_11 + we(nu)*norm*norm
231 8809 : norm = dznrm2(lapw%nv(2),zMat%data_c(lapw%nv(1) + atoms%nlotot + 1:, nu),1) ! dznrm2 returns the 2-norm of the vector.
232 9505 : q0_22 = q0_22 + we(nu)*norm*norm
233 : ENDDO
234 : ENDIF
235 696 : q0_11 = q0_11/cell%omtil
236 696 : q0_22 = q0_22/cell%omtil
237 : ELSE
238 6626 : IF (zmat%l_real) THEN
239 43162 : DO nu = 1, ne
240 40172 : norm = dnrm2(lapw%nv(jspin),zMat%data_r(:, nu),1) ! dnrm2 returns the 2-norm of the vector.
241 43162 : q0 = q0 + we(nu)*norm*norm
242 : ENDDO
243 : ELSE
244 51456 : DO nu = 1, ne
245 47820 : norm = dznrm2(lapw%nv(jspin),zMat%data_c(:, nu),1) ! dznrm2 returns the 2-norm of the vector.
246 51456 : q0 = q0 + we(nu)*norm*norm
247 : ENDDO
248 : ENDIF
249 6626 : q0 = q0/cell%omtil
250 : ENDIF
251 :
252 7322 : IF ((noco%l_noco).AND.(ikpt.LE.fmpi%isize)) THEN
253 85894 : dos%qis = 0.0
254 : END IF
255 : END IF
256 :
257 104123 : wtf(:ne) = we(:ne)/cell%omtil
258 7322 : IF (l_dfpt) wtf1(:ne) = we1(:ne)/cell%omtil
259 :
260 : ! LOOP OVER OCCUPIED STATES
261 104123 : DO nu = 1, ne
262 :
263 : ! Do inverse FFT of state and add related charge density to overall charge density on real-space mesh.
264 104123 : IF (noco%l_noco) THEN
265 8809 : forw = .FALSE.
266 8809 : IF (noco%l_ss) THEN
267 136 : CALL state%putComplexStateOnGrid(lapw, 1, zMat%data_c(1:lapw%nv(1),nu))
268 136 : CALL stateB%putComplexStateOnGrid(lapw, 2, zMat%data_c(lapw%nv(1) + atoms%nlotot+1:lapw%nv(1) + atoms%nlotot+lapw%nv(2),nu))
269 136 : IF (l_dfpt) THEN
270 0 : CALL stateq%putComplexStateOnGrid(lapwq, 1, zMat1%data_c(1:lapwq%nv(1),nu))
271 0 : CALL stateBq%putComplexStateOnGrid(lapwq, 2, zMat1%data_c(lapwq%nv(1) + atoms%nlotot+1:lapwq%nv(1) + atoms%nlotot+lapwq%nv(2),nu))
272 : END IF
273 : ELSE
274 8673 : CALL state%putComplexStateOnGrid(lapw, jspin, zMat%data_c(1:lapw%nv(jspin),nu))
275 8673 : CALL stateB%putComplexStateOnGrid(lapw, jspin, zMat%data_c(lapw%nv(1) + atoms%nlotot+1:lapw%nv(1) + atoms%nlotot+lapw%nv(jspin),nu))
276 8673 : IF (l_dfpt) THEN
277 0 : CALL stateq%putComplexStateOnGrid(lapwq, jspin, zMat1%data_c(1:lapwq%nv(1),nu))
278 0 : CALL stateBq%putComplexStateOnGrid(lapwq, jspin, zMat1%data_c(lapwq%nv(1) + atoms%nlotot+1:lapwq%nv(1) + atoms%nlotot+lapwq%nv(jspin),nu))
279 : END IF
280 : END IF
281 :
282 8809 : CALL fft_interface(3, state%dimensions(:), state%grid, forw, stateIndices)
283 8809 : CALL fft_interface(3, stateB%dimensions(:), stateB%grid, forw, stateBIndices)
284 8809 : IF (l_dfpt) THEN
285 0 : CALL fft_interface(3, stateq%dimensions(:), stateq%grid, forw, stateqIndices)
286 0 : CALL fft_interface(3, stateBq%dimensions(:), stateBq%grid, forw, stateBqIndices)
287 : END IF
288 :
289 : ! In the non-collinear case rho becomes a hermitian 2x2
290 : ! matrix (rhomatGrid).
291 : IF (.NOT.l_dfpt) THEN
292 37617465 : DO ir = 0, rhomatGrid(1)%gridLength - 1
293 : !In this order: rho_11, rho_22, m_x/2, m_y/2
294 37608656 : rhomatGrid(1)%grid(ir) = rhomatGrid(1)%grid(ir) + wtf(nu) * ABS(state%grid(ir))**2
295 37608656 : rhomatGrid(2)%grid(ir) = rhomatGrid(2)%grid(ir) + wtf(nu) * ABS(stateB%grid(ir))**2
296 37608656 : rhomatGrid(3)%grid(ir) = rhomatGrid(3)%grid(ir) + wtf(nu) * (REAL(state%grid(ir))*REAL(stateB%grid(ir)) + AIMAG(state%grid(ir))*AIMAG(stateB%grid(ir)))
297 37617465 : rhomatGrid(4)%grid(ir) = rhomatGrid(4)%grid(ir) + wtf(nu) * (REAL(state%grid(ir))*AIMAG(stateB%grid(ir)) - AIMAG(state%grid(ir))*REAL(stateB%grid(ir)))
298 : END DO
299 : ELSE
300 : !TODO: This looks ultra different for DFPT.
301 : !TODO: Only touch this once the magic minus is fully consistent.
302 : END IF
303 :
304 : ! In a non-collinear calculation the interstitial charge
305 : ! cannot be calculated by a simple substraction if the
306 : ! muffin-tin (and vacuum) charge is know, because the
307 : ! total charge does not need to be one in each spin-
308 : ! channel. Thus it has to be calculated explicitly, if
309 : ! it is needed.
310 8809 : IF (banddos%dos .OR. banddos%vacdos .OR. input%cdinf.OR.banddos%band) THEN
311 1029024 : DO ir = 0, state%gridLength - 1
312 1029000 : state%grid(ir) = ABS(state%grid(ir))**2
313 1029024 : stateB%grid(ir) = ABS(stateB%grid(ir))**2
314 : END DO
315 :
316 24 : forw = .TRUE.
317 24 : CALL fft_interface(3, state%dimensions(:), state%grid, forw, fieldSphereIndices)
318 24 : CALL fft_interface(3, stateB%dimensions(:), stateB%grid, forw, fieldSphereIndices)
319 :
320 24 : CALL state%takeFieldFromGrid(stars, cwk, stateFFTRadius+0.0005)
321 236208 : DO istr = 1, stars%ng3
322 236208 : cwk(istr) = REAL(stars%nstr(istr))*cwk(istr)
323 : END DO
324 96384 : DO istr = 1, stars%ng3_fft
325 96360 : CALL pwint(stars, atoms, sym, cell, istr, x)
326 96360 : dos%qis(ev_list(nu), ikpt, 1) = dos%qis(ev_list(nu), ikpt, 1) + REAL(cwk(istr)*x)/cell%omtil
327 96384 : dos%qTot(ev_list(nu), ikpt, 1) = dos%qTot(ev_list(nu), ikpt, 1) + REAL(cwk(istr)*x)/cell%omtil
328 : ENDDO
329 :
330 24 : CALL stateB%takeFieldFromGrid(stars, cwk, stateFFTRadius+0.0005)
331 236208 : DO istr = 1, stars%ng3
332 236208 : cwk(istr) = REAL(stars%nstr(istr))*cwk(istr)
333 : END DO
334 96384 : DO istr = 1, stars%ng3_fft
335 96360 : CALL pwint(stars, atoms, sym, cell, istr, x)
336 96360 : dos%qis(ev_list(nu), ikpt, input%jspins) = dos%qis(ev_list(nu), ikpt, input%jspins) + REAL(cwk(istr)*x)/cell%omtil
337 105169 : dos%qTot(ev_list(nu), ikpt, input%jspins) = dos%qTot(ev_list(nu), ikpt, input%jspins) + REAL(cwk(istr)*x)/cell%omtil
338 : ENDDO
339 : ENDIF
340 :
341 : ELSE
342 87992 : CALL state%putStateOnGrid(lapw, jSpin, zMat, nu)
343 87992 : IF (l_dfpt) THEN
344 0 : CALL stateq%putStateOnGrid(lapwq, jspin, zMat1, nu)
345 0 : IF (l_minusq) THEN
346 0 : CALL statemq%putStateOnGrid(lapwmq, jspin, zMat1m, nu)
347 : END IF
348 : END IF
349 :
350 87992 : forw = .FALSE.
351 : ! The following FFT is a general complex to complex FFT
352 : ! For zmat%l_real this should be turned into a real to real FFT at some point
353 : ! Note: For the moment also no zero-indices for SpFFT provided
354 : !$acc data copy(state%grid)
355 87992 : CALL fft_interface(3, state%dimensions(:), state%grid, forw, stateIndices)
356 : !$acc end data
357 87992 : IF (l_dfpt) THEN
358 0 : CALL fft_interface(3, stateq%dimensions(:), stateq%grid, forw, stateqIndices)
359 0 : IF (l_minusq) THEN
360 0 : CALL fft_interface(3, statemq%dimensions(:), statemq%grid, forw, statemqIndices)
361 : END IF
362 : END IF
363 : IF (.NOT.l_dfpt) THEN
364 302641942 : DO ir = 0, chargeDen%gridLength - 1
365 302641942 : chargeDen%grid(ir) = chargeDen%grid(ir) + wtf(nu) * ABS(state%grid(ir))**2
366 : END DO
367 : ELSE
368 0 : IF (.NOT.l_minusq) THEN
369 0 : DO ir = 0, chargeDen%gridLength - 1
370 0 : chargeDen%grid(ir) = chargeDen%grid(ir) + wtf(nu) * 2 * CONJG(state%grid(ir)) * stateq%grid(ir)
371 : !chargeDen%grid(ir) = chargeDen%grid(ir) + wtf(nu) * 2 * CONJG(CONJG(state%grid(ir)) * stateq%grid(ir))
372 0 : IF (norm2(q_dfpt)<1e-8) chargeDen%grid(ir) = chargeDen%grid(ir) + wtf1(nu) * ABS(state%grid(ir))**2
373 : END DO
374 : ELSE
375 0 : DO ir = 0, chargeDen%gridLength - 1
376 : chargeDen%grid(ir) = chargeDen%grid(ir) + wtf(nu) * (CONJG(statemq%grid(ir)) * state%grid(ir) &
377 0 : + CONJG(state%grid(ir)) * stateq%grid(ir))
378 : !chargeDen%grid(ir) = chargeDen%grid(ir) + wtf(nu) * 2 * CONJG(CONJG(state%grid(ir)) * stateq%grid(ir))
379 0 : IF (norm2(q_dfpt)<1e-8) chargeDen%grid(ir) = chargeDen%grid(ir) + wtf1(nu) * ABS(state%grid(ir))**2
380 : END DO
381 : END IF
382 : END IF
383 :
384 87992 : IF (input%l_f) THEN
385 6998766 : DO ir = 0, ekinGrid%gridLength - 1
386 6998766 : ekinGrid%grid(ir) = ekinGrid%grid(ir) - wtf(nu) * eig(nu) * ABS(state%grid(ir))**2
387 : END DO
388 :
389 1568 : DO j = 1, 3
390 917808 : DO iv = 1, lapw%nv(jspin)
391 3666528 : xk = lapw%gvec(:, iv, jspin) + lapw%bkpt
392 : s = 0.0
393 3666528 : DO i = 1, 3
394 3666528 : s = s + xk(i)*cell%bmat(i, j)
395 : ENDDO
396 917808 : IF (zmat%l_real) THEN
397 747558 : tempState(iv) = s*zMat%data_r(iv, nu)
398 : ELSE
399 169074 : tempState(iv) = s*zMat%data_c(iv, nu)
400 : END IF
401 : END DO
402 1176 : CALL stateDeriv%putComplexStateOnGrid(lapw, jSpin, tempState)
403 1176 : CALL fft_interface(3, stateDeriv%dimensions(:), stateDeriv%grid, forw, stateIndices)
404 20996690 : DO ir = 0, ekinGrid%gridLength - 1
405 20996298 : ekinGrid%grid(ir) = ekinGrid%grid(ir) + wtf(nu) * 0.5 * ABS(stateDeriv%grid(ir))**2
406 : END DO
407 : END DO
408 : END IF
409 :
410 87992 : IF (noco%l_soc.AND.input%jspins.EQ.2) THEN
411 17556 : IF (banddos%dos .OR. banddos%vacdos .OR. input%cdinf.OR.banddos%band) THEN
412 0 : DO ir = 0, state%gridLength - 1
413 0 : state%grid(ir) = conjg(state%grid(ir)) * stepFct%grid(ir) * state%grid(ir)
414 : END DO
415 :
416 0 : forw = .TRUE.
417 0 : CALL fft_interface(3, state%dimensions(:), state%grid, forw, fieldSphereIndices)
418 :
419 0 : cwk = CMPLX(0.0,0.0)
420 0 : CALL state%takeFieldFromGrid(stars, cwk, stateFFTRadius+0.0005)
421 0 : starCharges = CMPLX(0.0,0.0)
422 0 : CALL pwint_all(stars,atoms,sym ,cell,1,stars%ng3,starCharges)
423 0 : starCharges(:) = starCharges(:) * cwk(:) * stars%nstr(:) / cell%omtil
424 0 : dos%qis(ev_list(nu), ikpt, jSpin) = dos%qis(ev_list(nu), ikpt, jSpin) + REAL(SUM(starCharges(:)))
425 0 : dos%qTot(ev_list(nu), ikpt, jSpin) = dos%qTot(ev_list(nu), ikpt, jSpin) + REAL(SUM(starCharges(:)))
426 : END IF
427 : END IF
428 : END IF
429 : END DO
430 :
431 : ! END OUTER LOOP OVER STATES NU
432 :
433 : ! Perform FFT transform of charge density to reciprocal space.
434 :
435 : ! In a collinear calculation pwden is called once per spin.
436 : ! However in a non-collinear calculation pwden is only called once
437 : ! and all four components of the density matrix (n_11 n_22 n_12
438 : ! n_21) have to be calculated at once.
439 7322 : ndens = 1
440 7322 : IF (noco%l_noco) ndens = 4
441 16732 : DO idens = 1, ndens
442 9410 : forw = .TRUE.
443 9410 : IF (noco%l_noco) THEN
444 2784 : CALL fft_interface(3, rhomatGrid(idens)%dimensions(:), rhomatGrid(idens)%grid, forw, fieldSphereIndices)
445 : ELSE
446 : ! The following FFT is a general complex to complex FFT
447 : ! For zmat%l_real this should be turned into a real to real FFT at some point
448 : ! Note: For the moment also no zero-indices for SpFFT provided
449 6626 : CALL fft_interface(3, chargeDen%dimensions(:), chargeDen%grid, forw, fieldSphereIndices)
450 6626 : IF (input%l_f) THEN
451 62 : CALL fft_interface(3, ekinGrid%dimensions(:), ekinGrid%grid, forw, fieldSphereIndices)
452 : END IF
453 : ENDIF
454 :
455 : ! collect stars from the fft-grid
456 20584172 : cwk = 0.0
457 20584172 : ecwk = 0.0
458 9410 : IF (noco%l_noco) THEN
459 2784 : CALL rhomatGrid(idens)%takeFieldFromGrid(stars, cwk, stateFFTRadius+0.0005)
460 : ELSE
461 : ! TODO: Shouldn't there be a starsq here for DFPT?
462 6626 : CALL chargeDen%takeFieldFromGrid(stars, cwk, stateFFTRadius+0.0005)
463 6626 : IF (input%l_f) THEN
464 62 : CALL ekinGrid%takeFieldFromGrid(stars, ecwk, stateFFTRadius+0.0005)
465 : END IF
466 : ENDIF
467 :
468 9410 : IF (input%l_useapw) THEN
469 0 : IF (input%l_f) THEN
470 0 : CALL force_b8(atoms, ecwk, stars, sym, cell, jspin, results%force, f_b8)
471 : ENDIF
472 : ENDIF
473 :
474 9410 : IF (.NOT.l_dfpt) THEN
475 : ! check charge neutralilty
476 9410 : IF ((idens .EQ. 1) .OR. (idens .EQ. 2)) THEN
477 8018 : IF (noco%l_noco) THEN
478 1392 : IF (idens .EQ. 1) THEN
479 696 : q0 = q0_11
480 : ELSE
481 696 : q0 = q0_22
482 : ENDIF
483 : ENDIF
484 8018 : IF (ABS(q0) .GT. 1.0e-9) THEN
485 7936 : IF (ABS(q0 - REAL(cwk(1)))/q0 .GT. tol_3) THEN
486 0 : WRITE (99, *) "XX:", ne, lapw%nv
487 0 : IF (zmat%l_real) THEN
488 0 : DO istr = 1, SIZE(zMat%data_r, 2)
489 0 : WRITE (99, *) "X:", istr, zMat%data_r(:, istr)
490 : ENDDO
491 : ELSE
492 0 : DO istr = 1, SIZE(zMat%data_c, 2)
493 0 : WRITE (99, *) "X:", istr, zMat%data_c(:, istr)
494 : ENDDO
495 : ENDIF
496 0 : WRITE (oUnit, '(''bad quality of charge density'',2f13.8)') q0, REAL(cwk(1))
497 0 : CALL juDFT_warn('pwden: bad quality of charge')
498 : ENDIF
499 : ENDIF
500 : ENDIF
501 : END IF
502 :
503 : ! add charge density to existing one
504 16732 : IF (idens .LE. 2) THEN
505 : ! add to spin-up or -down density (collinear & non-collinear)
506 8018 : ispin = jspin
507 8018 : IF (noco%l_noco) ispin = idens
508 : ! TODO: Shouldn't there be a starsq here for DFPT?
509 3791476 : DO istr = 1, stars%ng3_fft
510 3791476 : den%pw(istr, ispin) = den%pw(istr, ispin) + cwk(istr)
511 : ENDDO
512 1392 : ELSE IF (idens .EQ. 3) THEN
513 : ! add to off-diag. part of density matrix (only non-collinear)
514 1040324 : DO istr = 1, stars%ng3_fft
515 1040324 : den%pw(istr, 3) = den%pw(istr, 3) + cwk(istr)
516 : ENDDO
517 696 : IF (l_dfpt) THEN
518 : ! add to other off-diag. part of density matrix (only non-collinear)
519 0 : DO istr = 1, stars%ng3_fft
520 0 : den%pw(istr, 4) = den%pw(istr, 4) + cwk(istr)
521 : ENDDO
522 : END IF
523 : ELSE
524 : ! add to off-diag. part of density matrix (only non-collinear)
525 1040324 : DO istr = 1, stars%ng3_fft
526 1040324 : den%pw(istr, 3) = den%pw(istr, 3) - ImagUnit*cwk(istr)
527 : ! TODO: This is a magic minus. It should be + ImagUnit*cwk(istr)
528 : ENDDO
529 696 : IF (l_dfpt) THEN
530 : ! TODO: Only touch this once the magic minus is fully consistent.
531 0 : DO istr = 1, stars%ng3_fft
532 0 : den%pw(istr, 4) = den%pw(istr, 4) + ImagUnit*cwk(istr)
533 : ENDDO
534 : END IF
535 : ENDIF
536 : ENDDO
537 :
538 7322 : DEALLOCATE (cwk, ecwk)
539 :
540 7322 : CALL timestop("pwden")
541 :
542 7322 : END SUBROUTINE pwden
543 : END MODULE m_pwden
|
