Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 : ! **************************************************************************************************
8 : MODULE qs_rho0_types
9 :
10 : USE cp_units, ONLY: cp_unit_from_cp2k
11 : USE kinds, ONLY: dp
12 : USE mathconstants, ONLY: fourpi,&
13 : pi,&
14 : rootpi
15 : USE memory_utilities, ONLY: reallocate
16 : USE pw_types, ONLY: pw_c1d_gs_type,&
17 : pw_r3d_rs_type
18 : USE qs_grid_atom, ONLY: grid_atom_type
19 : USE qs_rho_atom_types, ONLY: rho_atom_coeff
20 : USE whittaker, ONLY: whittaker_c0a,&
21 : whittaker_ci
22 : #include "./base/base_uses.f90"
23 :
24 : IMPLICIT NONE
25 :
26 : PRIVATE
27 :
28 : ! *** Global parameters (only in this module)
29 :
30 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_rho0_types'
31 :
32 : ! *** Define multipole type ***
33 :
34 : ! **************************************************************************************************
35 : TYPE mpole_rho_atom
36 : REAL(dp), DIMENSION(:), POINTER :: Qlm_h => NULL(), &
37 : Qlm_s => NULL(), &
38 : Qlm_tot => NULL(), &
39 : Qlm_car => NULL()
40 : REAL(dp) :: Qlm_z = -1.0_dp
41 : REAL(dp), DIMENSION(2) :: Q0 = -1.0_dp
42 : END TYPE mpole_rho_atom
43 :
44 : ! **************************************************************************************************
45 : TYPE mpole_gau_overlap
46 : REAL(dp), DIMENSION(:, :, :), POINTER :: Qlm_gg => NULL()
47 : REAL(dp), DIMENSION(:, :), POINTER :: g0_h => NULL(), vg0_h => NULL()
48 : REAL(dp) :: rpgf0_h = -1.0_dp, rpgf0_s = -1.0_dp
49 : END TYPE mpole_gau_overlap
50 :
51 : ! **************************************************************************************************
52 : TYPE rho0_mpole_type
53 : TYPE(mpole_rho_atom), DIMENSION(:), POINTER :: mp_rho => NULL()
54 : TYPE(mpole_gau_overlap), DIMENSION(:), &
55 : POINTER :: mp_gau => NULL()
56 : REAL(dp) :: zet0_h = -1.0_dp, &
57 : total_rho0_h = -1.0_dp
58 : REAL(dp) :: max_rpgf0_s = -1.0_dp
59 : REAL(dp), DIMENSION(:), POINTER :: norm_g0l_h => NULL()
60 : INTEGER, DIMENSION(:), POINTER :: lmax0_kind => NULL()
61 : INTEGER :: lmax_0 = -1, igrid_zet0_s = -1
62 : TYPE(pw_r3d_rs_type), POINTER :: rho0_s_rs => NULL()
63 : TYPE(pw_c1d_gs_type), POINTER :: rho0_s_gs => NULL()
64 : END TYPE rho0_mpole_type
65 :
66 : ! **************************************************************************************************
67 : TYPE rho0_atom_type
68 : TYPE(rho_atom_coeff), POINTER :: rho0_rad_h => NULL(), &
69 : vrho0_rad_h => NULL()
70 : END TYPE rho0_atom_type
71 :
72 : ! Public Types
73 :
74 : PUBLIC :: mpole_rho_atom, mpole_gau_overlap, &
75 : rho0_atom_type, rho0_mpole_type
76 :
77 : ! Public Subroutine
78 :
79 : PUBLIC :: allocate_multipoles, allocate_rho0_mpole, &
80 : allocate_rho0_atom, allocate_rho0_atom_rad, &
81 : deallocate_rho0_atom, deallocate_rho0_mpole, &
82 : calculate_g0, get_rho0_mpole, initialize_mpole_rho, &
83 : write_rho0_info
84 :
85 : CONTAINS
86 :
87 : ! **************************************************************************************************
88 : !> \brief ...
89 : !> \param mp_rho ...
90 : !> \param natom ...
91 : !> \param mp_gau ...
92 : !> \param nkind ...
93 : ! **************************************************************************************************
94 1792 : SUBROUTINE allocate_multipoles(mp_rho, natom, mp_gau, nkind)
95 :
96 : TYPE(mpole_rho_atom), DIMENSION(:), POINTER :: mp_rho
97 : INTEGER, INTENT(IN) :: natom
98 : TYPE(mpole_gau_overlap), DIMENSION(:), POINTER :: mp_gau
99 : INTEGER, INTENT(IN) :: nkind
100 :
101 : INTEGER :: iat, ikind
102 :
103 1792 : IF (ASSOCIATED(mp_rho)) THEN
104 0 : CALL deallocate_mpole_rho(mp_rho)
105 : END IF
106 :
107 15092 : ALLOCATE (mp_rho(natom))
108 :
109 7924 : DO iat = 1, natom
110 6132 : NULLIFY (mp_rho(iat)%Qlm_h)
111 6132 : NULLIFY (mp_rho(iat)%Qlm_s)
112 6132 : NULLIFY (mp_rho(iat)%Qlm_tot)
113 7924 : NULLIFY (mp_rho(iat)%Qlm_car)
114 : END DO
115 :
116 1792 : IF (ASSOCIATED(mp_gau)) THEN
117 0 : CALL deallocate_mpole_gau(mp_gau)
118 : END IF
119 :
120 9122 : ALLOCATE (mp_gau(nkind))
121 :
122 5538 : DO ikind = 1, nkind
123 3746 : NULLIFY (mp_gau(ikind)%Qlm_gg)
124 3746 : NULLIFY (mp_gau(ikind)%g0_h)
125 3746 : NULLIFY (mp_gau(ikind)%vg0_h)
126 3746 : mp_gau(ikind)%rpgf0_h = 0.0_dp
127 5538 : mp_gau(ikind)%rpgf0_s = 0.0_dp
128 : END DO
129 :
130 1792 : END SUBROUTINE allocate_multipoles
131 :
132 : ! **************************************************************************************************
133 : !> \brief ...
134 : !> \param rho0_set ...
135 : !> \param natom ...
136 : ! **************************************************************************************************
137 1792 : SUBROUTINE allocate_rho0_atom(rho0_set, natom)
138 :
139 : TYPE(rho0_atom_type), DIMENSION(:), POINTER :: rho0_set
140 : INTEGER, INTENT(IN) :: natom
141 :
142 : INTEGER :: iat
143 :
144 1792 : IF (ASSOCIATED(rho0_set)) THEN
145 0 : CALL deallocate_rho0_atom(rho0_set)
146 : END IF
147 :
148 11508 : ALLOCATE (rho0_set(natom))
149 :
150 7924 : DO iat = 1, natom
151 6132 : NULLIFY (rho0_set(iat)%rho0_rad_h)
152 7924 : NULLIFY (rho0_set(iat)%vrho0_rad_h)
153 : END DO
154 :
155 1792 : END SUBROUTINE allocate_rho0_atom
156 :
157 : ! **************************************************************************************************
158 : !> \brief ...
159 : !> \param rho0_atom ...
160 : !> \param nr ...
161 : !> \param nchannels ...
162 : ! **************************************************************************************************
163 6132 : SUBROUTINE allocate_rho0_atom_rad(rho0_atom, nr, nchannels)
164 :
165 : TYPE(rho0_atom_type), INTENT(OUT) :: rho0_atom
166 : INTEGER, INTENT(IN) :: nr, nchannels
167 :
168 6132 : ALLOCATE (rho0_atom%rho0_rad_h)
169 :
170 : NULLIFY (rho0_atom%rho0_rad_h%r_coef)
171 24528 : ALLOCATE (rho0_atom%rho0_rad_h%r_coef(1:nr, 1:nchannels))
172 2394128 : rho0_atom%rho0_rad_h%r_coef = 0.0_dp
173 :
174 6132 : ALLOCATE (rho0_atom%vrho0_rad_h)
175 :
176 : NULLIFY (rho0_atom%vrho0_rad_h%r_coef)
177 24528 : ALLOCATE (rho0_atom%vrho0_rad_h%r_coef(1:nr, 1:nchannels))
178 2394128 : rho0_atom%vrho0_rad_h%r_coef = 0.0_dp
179 :
180 6132 : END SUBROUTINE allocate_rho0_atom_rad
181 :
182 : ! **************************************************************************************************
183 : !> \brief ...
184 : !> \param rho0 ...
185 : ! **************************************************************************************************
186 1792 : SUBROUTINE allocate_rho0_mpole(rho0)
187 :
188 : TYPE(rho0_mpole_type), POINTER :: rho0
189 :
190 1792 : IF (ASSOCIATED(rho0)) THEN
191 0 : CALL deallocate_rho0_mpole(rho0)
192 : END IF
193 :
194 1792 : ALLOCATE (rho0)
195 :
196 : NULLIFY (rho0%mp_rho)
197 : NULLIFY (rho0%mp_gau)
198 : NULLIFY (rho0%norm_g0l_h)
199 : NULLIFY (rho0%lmax0_kind)
200 : NULLIFY (rho0%rho0_s_rs)
201 : NULLIFY (rho0%rho0_s_gs)
202 :
203 1792 : END SUBROUTINE allocate_rho0_mpole
204 :
205 : ! **************************************************************************************************
206 : !> \brief ...
207 : !> \param rho0_mpole ...
208 : !> \param grid_atom ...
209 : !> \param ik ...
210 : ! **************************************************************************************************
211 3746 : SUBROUTINE calculate_g0(rho0_mpole, grid_atom, ik)
212 :
213 : TYPE(rho0_mpole_type), POINTER :: rho0_mpole
214 : TYPE(grid_atom_type), POINTER :: grid_atom
215 : INTEGER, INTENT(IN) :: ik
216 :
217 : INTEGER :: ir, l, lmax, nr
218 : REAL(dp) :: c1, prefactor, root_z_h, z_h
219 3746 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: erf_z_h, gexp, gh_tmp, int1, int2
220 :
221 3746 : nr = grid_atom%nr
222 3746 : lmax = rho0_mpole%lmax0_kind(ik)
223 3746 : z_h = rho0_mpole%zet0_h
224 3746 : root_z_h = SQRT(z_h)
225 :
226 : ! Allocate g0
227 3746 : CALL reallocate(rho0_mpole%mp_gau(ik)%g0_h, 1, nr, 0, lmax)
228 3746 : CALL reallocate(rho0_mpole%mp_gau(ik)%vg0_h, 1, nr, 0, lmax)
229 :
230 26222 : ALLOCATE (gexp(nr), gh_tmp(nr), erf_z_h(nr), int1(nr), int2(nr))
231 :
232 195966 : gh_tmp(1:nr) = EXP(-z_h*grid_atom%rad2(1:nr))
233 :
234 195966 : DO ir = 1, nr
235 195966 : erf_z_h(ir) = erf(grid_atom%rad(ir)*root_z_h)
236 : END DO
237 :
238 195966 : DO ir = 1, nr
239 195966 : IF (gh_tmp(ir) < 1.0E-30_dp) gh_tmp(ir) = 0.0_dp
240 : END DO
241 :
242 195966 : gexp(1:nr) = gh_tmp(1:nr)
243 : rho0_mpole%mp_gau(ik)%g0_h(1:nr, 0) = gh_tmp(1:nr)* &
244 195966 : rho0_mpole%norm_g0l_h(0)
245 3746 : CALL whittaker_c0a(int1, grid_atom%rad, gh_tmp, erf_z_h, z_h, 0, 0, nr)
246 3746 : CALL whittaker_ci(int2, grid_atom%rad, gh_tmp, z_h, 0, nr)
247 :
248 3746 : prefactor = fourpi*rho0_mpole%norm_g0l_h(0)
249 :
250 3746 : c1 = SQRT(pi*pi*pi/(z_h*z_h*z_h))*rho0_mpole%norm_g0l_h(0)
251 :
252 195966 : DO ir = 1, nr
253 195966 : rho0_mpole%mp_gau(ik)%vg0_h(ir, 0) = c1*erf_z_h(ir)*grid_atom%oorad2l(ir, 1)
254 : END DO
255 :
256 10234 : DO l = 1, lmax
257 339848 : gh_tmp(1:nr) = gh_tmp(1:nr)*grid_atom%rad(1:nr)
258 : rho0_mpole%mp_gau(ik)%g0_h(1:nr, l) = gh_tmp(1:nr)* &
259 339848 : rho0_mpole%norm_g0l_h(l)
260 :
261 6488 : prefactor = fourpi/(2.0_dp*l + 1.0_dp)*rho0_mpole%norm_g0l_h(l)
262 6488 : CALL whittaker_c0a(int1, grid_atom%rad, gexp, erf_z_h, z_h, l, l, nr)
263 343594 : DO ir = 1, nr
264 : rho0_mpole%mp_gau(ik)%vg0_h(ir, l) = prefactor*(int1(ir) + &
265 339848 : int2(ir)*grid_atom%rad2l(ir, l))
266 : END DO
267 :
268 : END DO ! l
269 :
270 3746 : DEALLOCATE (gexp, erf_z_h, gh_tmp, int1, int2)
271 3746 : END SUBROUTINE calculate_g0
272 :
273 : ! **************************************************************************************************
274 : !> \brief ...
275 : !> \param mp_gau ...
276 : ! **************************************************************************************************
277 1792 : SUBROUTINE deallocate_mpole_gau(mp_gau)
278 :
279 : TYPE(mpole_gau_overlap), DIMENSION(:), POINTER :: mp_gau
280 :
281 : INTEGER :: ikind, nkind
282 :
283 1792 : nkind = SIZE(mp_gau)
284 :
285 5538 : DO ikind = 1, nkind
286 :
287 3746 : IF (ASSOCIATED(mp_gau(ikind)%Qlm_gg)) THEN
288 3402 : DEALLOCATE (mp_gau(ikind)%Qlm_gg)
289 : END IF
290 :
291 3746 : DEALLOCATE (mp_gau(ikind)%g0_h)
292 :
293 5538 : DEALLOCATE (mp_gau(ikind)%vg0_h)
294 : END DO
295 :
296 1792 : DEALLOCATE (mp_gau)
297 :
298 1792 : END SUBROUTINE deallocate_mpole_gau
299 :
300 : ! **************************************************************************************************
301 : !> \brief ...
302 : !> \param mp_rho ...
303 : ! **************************************************************************************************
304 1792 : SUBROUTINE deallocate_mpole_rho(mp_rho)
305 :
306 : TYPE(mpole_rho_atom), DIMENSION(:), POINTER :: mp_rho
307 :
308 : INTEGER :: iat, natom
309 :
310 1792 : natom = SIZE(mp_rho)
311 :
312 7924 : DO iat = 1, natom
313 6132 : DEALLOCATE (mp_rho(iat)%Qlm_h)
314 6132 : DEALLOCATE (mp_rho(iat)%Qlm_s)
315 6132 : DEALLOCATE (mp_rho(iat)%Qlm_tot)
316 7924 : DEALLOCATE (mp_rho(iat)%Qlm_car)
317 : END DO
318 :
319 1792 : DEALLOCATE (mp_rho)
320 :
321 1792 : END SUBROUTINE deallocate_mpole_rho
322 :
323 : ! **************************************************************************************************
324 : !> \brief ...
325 : !> \param rho0_atom_set ...
326 : ! **************************************************************************************************
327 1792 : SUBROUTINE deallocate_rho0_atom(rho0_atom_set)
328 :
329 : TYPE(rho0_atom_type), DIMENSION(:), POINTER :: rho0_atom_set
330 :
331 : INTEGER :: iat, natom
332 :
333 1792 : IF (ASSOCIATED(rho0_atom_set)) THEN
334 :
335 1792 : natom = SIZE(rho0_atom_set)
336 :
337 7924 : DO iat = 1, natom
338 6132 : IF (ASSOCIATED(rho0_atom_set(iat)%rho0_rad_h)) THEN
339 6132 : DEALLOCATE (rho0_atom_set(iat)%rho0_rad_h%r_coef)
340 6132 : DEALLOCATE (rho0_atom_set(iat)%rho0_rad_h)
341 : END IF
342 7924 : IF (ASSOCIATED(rho0_atom_set(iat)%vrho0_rad_h)) THEN
343 6132 : DEALLOCATE (rho0_atom_set(iat)%vrho0_rad_h%r_coef)
344 6132 : DEALLOCATE (rho0_atom_set(iat)%vrho0_rad_h)
345 : END IF
346 : END DO
347 :
348 1792 : DEALLOCATE (rho0_atom_set)
349 : ELSE
350 : CALL cp_abort(__LOCATION__, &
351 : "The pointer rho0_atom_set is not associated and "// &
352 0 : "cannot be deallocated")
353 : END IF
354 :
355 1792 : END SUBROUTINE deallocate_rho0_atom
356 : ! **************************************************************************************************
357 : !> \brief ...
358 : !> \param rho0 ...
359 : ! **************************************************************************************************
360 1792 : SUBROUTINE deallocate_rho0_mpole(rho0)
361 :
362 : TYPE(rho0_mpole_type), POINTER :: rho0
363 :
364 1792 : IF (ASSOCIATED(rho0)) THEN
365 :
366 1792 : IF (ASSOCIATED(rho0%mp_gau)) CALL deallocate_mpole_gau(rho0%mp_gau)
367 :
368 1792 : IF (ASSOCIATED(rho0%mp_rho)) CALL deallocate_mpole_rho(rho0%mp_rho)
369 :
370 1792 : IF (ASSOCIATED(rho0%lmax0_kind)) THEN
371 1792 : DEALLOCATE (rho0%lmax0_kind)
372 : END IF
373 :
374 1792 : IF (ASSOCIATED(rho0%norm_g0l_h)) THEN
375 1792 : DEALLOCATE (rho0%norm_g0l_h)
376 : END IF
377 :
378 1792 : IF (ASSOCIATED(rho0%rho0_s_rs)) THEN
379 1792 : CALL rho0%rho0_s_rs%release()
380 1792 : DEALLOCATE (rho0%rho0_s_rs)
381 : END IF
382 :
383 1792 : IF (ASSOCIATED(rho0%rho0_s_gs)) THEN
384 1792 : CALL rho0%rho0_s_gs%release()
385 1792 : DEALLOCATE (rho0%rho0_s_gs)
386 :
387 : END IF
388 1792 : DEALLOCATE (rho0)
389 : ELSE
390 : CALL cp_abort(__LOCATION__, &
391 : "The pointer rho0 is not associated and "// &
392 0 : "cannot be deallocated")
393 : END IF
394 :
395 1792 : END SUBROUTINE deallocate_rho0_mpole
396 :
397 : ! **************************************************************************************************
398 : !> \brief ...
399 : !> \param rho0_mpole ...
400 : !> \param g0_h ...
401 : !> \param vg0_h ...
402 : !> \param iat ...
403 : !> \param ikind ...
404 : !> \param lmax_0 ...
405 : !> \param l0_ikind ...
406 : !> \param mp_gau_ikind ...
407 : !> \param mp_rho ...
408 : !> \param norm_g0l_h ...
409 : !> \param Qlm_gg ...
410 : !> \param Qlm_car ...
411 : !> \param Qlm_tot ...
412 : !> \param zet0_h ...
413 : !> \param igrid_zet0_s ...
414 : !> \param rpgf0_h ...
415 : !> \param rpgf0_s ...
416 : !> \param max_rpgf0_s ...
417 : !> \param rho0_s_rs ...
418 : !> \param rho0_s_gs ...
419 : ! **************************************************************************************************
420 219277 : SUBROUTINE get_rho0_mpole(rho0_mpole, g0_h, vg0_h, iat, ikind, lmax_0, l0_ikind, &
421 : mp_gau_ikind, mp_rho, norm_g0l_h, &
422 : Qlm_gg, Qlm_car, Qlm_tot, &
423 : zet0_h, igrid_zet0_s, rpgf0_h, rpgf0_s, &
424 : max_rpgf0_s, rho0_s_rs, rho0_s_gs)
425 :
426 : TYPE(rho0_mpole_type), POINTER :: rho0_mpole
427 : REAL(dp), DIMENSION(:, :), OPTIONAL, POINTER :: g0_h, vg0_h
428 : INTEGER, INTENT(IN), OPTIONAL :: iat, ikind
429 : INTEGER, INTENT(OUT), OPTIONAL :: lmax_0, l0_ikind
430 : TYPE(mpole_gau_overlap), OPTIONAL, POINTER :: mp_gau_ikind
431 : TYPE(mpole_rho_atom), DIMENSION(:), OPTIONAL, &
432 : POINTER :: mp_rho
433 : REAL(dp), DIMENSION(:), OPTIONAL, POINTER :: norm_g0l_h
434 : REAL(dp), DIMENSION(:, :, :), OPTIONAL, POINTER :: Qlm_gg
435 : REAL(dp), DIMENSION(:), OPTIONAL, POINTER :: Qlm_car, Qlm_tot
436 : REAL(dp), INTENT(OUT), OPTIONAL :: zet0_h
437 : INTEGER, INTENT(OUT), OPTIONAL :: igrid_zet0_s
438 : REAL(dp), INTENT(OUT), OPTIONAL :: rpgf0_h, rpgf0_s, max_rpgf0_s
439 : TYPE(pw_r3d_rs_type), OPTIONAL, POINTER :: rho0_s_rs
440 : TYPE(pw_c1d_gs_type), OPTIONAL, POINTER :: rho0_s_gs
441 :
442 219277 : IF (ASSOCIATED(rho0_mpole)) THEN
443 :
444 219277 : IF (PRESENT(lmax_0)) lmax_0 = rho0_mpole%lmax_0
445 219277 : IF (PRESENT(mp_rho)) mp_rho => rho0_mpole%mp_rho
446 219277 : IF (PRESENT(norm_g0l_h)) norm_g0l_h => rho0_mpole%norm_g0l_h
447 219277 : IF (PRESENT(zet0_h)) zet0_h = rho0_mpole%zet0_h
448 219277 : IF (PRESENT(igrid_zet0_s)) igrid_zet0_s = rho0_mpole%igrid_zet0_s
449 219277 : IF (PRESENT(max_rpgf0_s)) max_rpgf0_s = rho0_mpole%max_rpgf0_s
450 219277 : IF (PRESENT(rho0_s_rs)) rho0_s_rs => rho0_mpole%rho0_s_rs
451 219277 : IF (PRESENT(rho0_s_gs)) rho0_s_gs => rho0_mpole%rho0_s_gs
452 :
453 219277 : IF (PRESENT(ikind)) THEN
454 124444 : IF (PRESENT(l0_ikind)) l0_ikind = rho0_mpole%lmax0_kind(ikind)
455 124444 : IF (PRESENT(mp_gau_ikind)) mp_gau_ikind => rho0_mpole%mp_gau(ikind)
456 124444 : IF (PRESENT(g0_h)) g0_h => rho0_mpole%mp_gau(ikind)%g0_h
457 124444 : IF (PRESENT(vg0_h)) vg0_h => rho0_mpole%mp_gau(ikind)%vg0_h
458 124444 : IF (PRESENT(Qlm_gg)) Qlm_gg => rho0_mpole%mp_gau(ikind)%Qlm_gg
459 124444 : IF (PRESENT(rpgf0_h)) rpgf0_h = rho0_mpole%mp_gau(ikind)%rpgf0_h
460 124444 : IF (PRESENT(rpgf0_s)) rpgf0_s = rho0_mpole%mp_gau(ikind)%rpgf0_s
461 : END IF
462 219277 : IF (PRESENT(iat)) THEN
463 47251 : IF (PRESENT(Qlm_car)) Qlm_car => rho0_mpole%mp_rho(iat)%Qlm_car
464 47251 : IF (PRESENT(Qlm_tot)) Qlm_tot => rho0_mpole%mp_rho(iat)%Qlm_tot
465 : END IF
466 :
467 : ELSE
468 0 : CPABORT("The pointer rho0_mpole is not associated")
469 : END IF
470 :
471 219277 : END SUBROUTINE get_rho0_mpole
472 :
473 : ! **************************************************************************************************
474 : !> \brief ...
475 : !> \param mp_rho ...
476 : !> \param nchan_s ...
477 : !> \param nchan_c ...
478 : !> \param zeff ...
479 : ! **************************************************************************************************
480 6132 : SUBROUTINE initialize_mpole_rho(mp_rho, nchan_s, nchan_c, zeff)
481 :
482 : TYPE(mpole_rho_atom) :: mp_rho
483 : INTEGER, INTENT(IN) :: nchan_s, nchan_c
484 : REAL(KIND=dp), INTENT(IN) :: zeff
485 :
486 6132 : CALL reallocate(mp_rho%Qlm_h, 1, nchan_s)
487 6132 : CALL reallocate(mp_rho%Qlm_s, 1, nchan_s)
488 6132 : CALL reallocate(mp_rho%Qlm_tot, 1, nchan_s)
489 6132 : CALL reallocate(mp_rho%Qlm_car, 1, nchan_c)
490 :
491 52088 : mp_rho%Qlm_h = 0.0_dp
492 52088 : mp_rho%Qlm_s = 0.0_dp
493 52088 : mp_rho%Qlm_tot = 0.0_dp
494 57934 : mp_rho%Qlm_car = 0.0_dp
495 6132 : mp_rho%Qlm_z = -2.0_dp*rootpi*Zeff
496 18396 : mp_rho%Q0 = 0.0_dp
497 :
498 6132 : END SUBROUTINE initialize_mpole_rho
499 :
500 : ! **************************************************************************************************
501 : !> \brief ...
502 : !> \param rho0_mpole ...
503 : !> \param unit_str ...
504 : !> \param output_unit ...
505 : ! **************************************************************************************************
506 1 : SUBROUTINE write_rho0_info(rho0_mpole, unit_str, output_unit)
507 :
508 : TYPE(rho0_mpole_type), POINTER :: rho0_mpole
509 : CHARACTER(LEN=*), INTENT(IN) :: unit_str
510 : INTEGER, INTENT(in) :: output_unit
511 :
512 : INTEGER :: ikind, l, nkind
513 : REAL(dp) :: conv
514 :
515 1 : IF (ASSOCIATED(rho0_mpole)) THEN
516 1 : conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
517 :
518 : WRITE (UNIT=output_unit, FMT="(/,T5,A,/)") &
519 1 : "*** Compensation density charges data set ***"
520 : WRITE (UNIT=output_unit, FMT="(T2,A,T35,f16.10)") &
521 1 : "- Rho0 exponent :", rho0_mpole%zet0_h
522 : WRITE (UNIT=output_unit, FMT="(T2,A,T35,I5)") &
523 1 : "- Global max l :", rho0_mpole%lmax_0
524 :
525 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
526 1 : "- Normalization constants for g0"
527 2 : DO l = 0, rho0_mpole%lmax_0
528 : WRITE (UNIT=output_unit, FMT="(T20,A,T31,I2,T38,A,f15.5)") &
529 2 : "ang. mom.= ", l, " hard= ", rho0_mpole%norm_g0l_h(l)
530 : END DO
531 :
532 1 : nkind = SIZE(rho0_mpole%lmax0_kind, 1)
533 2 : DO ikind = 1, nkind
534 : WRITE (UNIT=output_unit, FMT="(/,T2,A,T55,I2)") &
535 : "- rho0 max L and radii in "//TRIM(unit_str)// &
536 1 : " for the atom kind :", ikind
537 :
538 : WRITE (UNIT=output_unit, FMT="(T2,T20,A,T55,I5)") &
539 1 : "=> l max :", rho0_mpole%lmax0_kind(ikind)
540 :
541 : WRITE (UNIT=output_unit, FMT="(T2,T20,A,T55,f20.10)") &
542 1 : "=> max radius of g0: ", &
543 3 : rho0_mpole%mp_gau(ikind)%rpgf0_h*conv
544 : END DO ! ikind
545 :
546 : ELSE
547 : WRITE (UNIT=output_unit, FMT="(/,T5,A,/)") &
548 0 : ' WARNING: I cannot print rho0, it is not associated'
549 : END IF
550 :
551 1 : END SUBROUTINE write_rho0_info
552 0 : END MODULE qs_rho0_types
|
