Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 :
65 : ! CNEO nuclear charge density related stuff.
66 : ! These are put here because it is preferred that function rho0_s_grid_create
67 : ! can initialize PW for both rho0 and nuclear rho1s. This is akin to the behavior
68 : ! that init_rho0 also initializes data structures for CNEO nuclear charge densities,
69 : ! such that code modifications are mimimal for places where init_rho0 and
70 : ! rho0_s_grid_create are called.
71 : LOGICAL :: do_cneo = .FALSE. ! if CNEO potential is used
72 : REAL(dp) :: tot_rhoz_cneo_s = 0.0_dp ! soft nuclear charge on the local grid
73 : TYPE(pw_r3d_rs_type), POINTER :: rhoz_cneo_s_rs => Null() ! soft nuclear charge density, r-space
74 : TYPE(pw_c1d_gs_type), POINTER :: rhoz_cneo_s_gs => Null() ! soft nuclear charge density, g-space
75 : END TYPE rho0_mpole_type
76 :
77 : ! **************************************************************************************************
78 : TYPE rho0_atom_type
79 : TYPE(rho_atom_coeff), POINTER :: rho0_rad_h => NULL(), &
80 : vrho0_rad_h => NULL()
81 : END TYPE rho0_atom_type
82 :
83 : ! Public Types
84 :
85 : PUBLIC :: mpole_rho_atom, mpole_gau_overlap, &
86 : rho0_atom_type, rho0_mpole_type
87 :
88 : ! Public Subroutine
89 :
90 : PUBLIC :: allocate_multipoles, allocate_rho0_mpole, &
91 : allocate_rho0_atom, allocate_rho0_atom_rad, &
92 : deallocate_rho0_atom, deallocate_rho0_mpole, &
93 : calculate_g0, get_rho0_mpole, initialize_mpole_rho, &
94 : write_rho0_info
95 :
96 : CONTAINS
97 :
98 : ! **************************************************************************************************
99 : !> \brief ...
100 : !> \param mp_rho ...
101 : !> \param natom ...
102 : !> \param mp_gau ...
103 : !> \param nkind ...
104 : ! **************************************************************************************************
105 2168 : SUBROUTINE allocate_multipoles(mp_rho, natom, mp_gau, nkind)
106 :
107 : TYPE(mpole_rho_atom), DIMENSION(:), POINTER :: mp_rho
108 : INTEGER, INTENT(IN) :: natom
109 : TYPE(mpole_gau_overlap), DIMENSION(:), POINTER :: mp_gau
110 : INTEGER, INTENT(IN) :: nkind
111 :
112 : INTEGER :: iat, ikind
113 :
114 2168 : IF (ASSOCIATED(mp_rho)) THEN
115 0 : CALL deallocate_mpole_rho(mp_rho)
116 : END IF
117 :
118 17888 : ALLOCATE (mp_rho(natom))
119 :
120 9216 : DO iat = 1, natom
121 7048 : NULLIFY (mp_rho(iat)%Qlm_h)
122 7048 : NULLIFY (mp_rho(iat)%Qlm_s)
123 7048 : NULLIFY (mp_rho(iat)%Qlm_tot)
124 9216 : NULLIFY (mp_rho(iat)%Qlm_car)
125 : END DO
126 :
127 2168 : IF (ASSOCIATED(mp_gau)) THEN
128 0 : CALL deallocate_mpole_gau(mp_gau)
129 : END IF
130 :
131 10784 : ALLOCATE (mp_gau(nkind))
132 :
133 6448 : DO ikind = 1, nkind
134 4280 : NULLIFY (mp_gau(ikind)%Qlm_gg)
135 4280 : NULLIFY (mp_gau(ikind)%g0_h)
136 4280 : NULLIFY (mp_gau(ikind)%vg0_h)
137 4280 : mp_gau(ikind)%rpgf0_h = 0.0_dp
138 6448 : mp_gau(ikind)%rpgf0_s = 0.0_dp
139 : END DO
140 :
141 2168 : END SUBROUTINE allocate_multipoles
142 :
143 : ! **************************************************************************************************
144 : !> \brief ...
145 : !> \param rho0_set ...
146 : !> \param natom ...
147 : ! **************************************************************************************************
148 2168 : SUBROUTINE allocate_rho0_atom(rho0_set, natom)
149 :
150 : TYPE(rho0_atom_type), DIMENSION(:), POINTER :: rho0_set
151 : INTEGER, INTENT(IN) :: natom
152 :
153 : INTEGER :: iat
154 :
155 2168 : IF (ASSOCIATED(rho0_set)) THEN
156 0 : CALL deallocate_rho0_atom(rho0_set)
157 : END IF
158 :
159 13552 : ALLOCATE (rho0_set(natom))
160 :
161 9216 : DO iat = 1, natom
162 7048 : NULLIFY (rho0_set(iat)%rho0_rad_h)
163 9216 : NULLIFY (rho0_set(iat)%vrho0_rad_h)
164 : END DO
165 :
166 2168 : END SUBROUTINE allocate_rho0_atom
167 :
168 : ! **************************************************************************************************
169 : !> \brief ...
170 : !> \param rho0_atom ...
171 : !> \param nr ...
172 : !> \param nchannels ...
173 : ! **************************************************************************************************
174 7048 : SUBROUTINE allocate_rho0_atom_rad(rho0_atom, nr, nchannels)
175 :
176 : TYPE(rho0_atom_type), INTENT(OUT) :: rho0_atom
177 : INTEGER, INTENT(IN) :: nr, nchannels
178 :
179 7048 : ALLOCATE (rho0_atom%rho0_rad_h)
180 :
181 : NULLIFY (rho0_atom%rho0_rad_h%r_coef)
182 28192 : ALLOCATE (rho0_atom%rho0_rad_h%r_coef(1:nr, 1:nchannels))
183 2760816 : rho0_atom%rho0_rad_h%r_coef = 0.0_dp
184 :
185 7048 : ALLOCATE (rho0_atom%vrho0_rad_h)
186 :
187 : NULLIFY (rho0_atom%vrho0_rad_h%r_coef)
188 28192 : ALLOCATE (rho0_atom%vrho0_rad_h%r_coef(1:nr, 1:nchannels))
189 2760816 : rho0_atom%vrho0_rad_h%r_coef = 0.0_dp
190 :
191 7048 : END SUBROUTINE allocate_rho0_atom_rad
192 :
193 : ! **************************************************************************************************
194 : !> \brief ...
195 : !> \param rho0 ...
196 : ! **************************************************************************************************
197 2168 : SUBROUTINE allocate_rho0_mpole(rho0)
198 :
199 : TYPE(rho0_mpole_type), POINTER :: rho0
200 :
201 2168 : IF (ASSOCIATED(rho0)) THEN
202 0 : CALL deallocate_rho0_mpole(rho0)
203 : END IF
204 :
205 2168 : ALLOCATE (rho0)
206 :
207 : NULLIFY (rho0%mp_rho)
208 : NULLIFY (rho0%mp_gau)
209 : NULLIFY (rho0%norm_g0l_h)
210 : NULLIFY (rho0%lmax0_kind)
211 : NULLIFY (rho0%rho0_s_rs)
212 : NULLIFY (rho0%rho0_s_gs)
213 :
214 2168 : END SUBROUTINE allocate_rho0_mpole
215 :
216 : ! **************************************************************************************************
217 : !> \brief ...
218 : !> \param rho0_mpole ...
219 : !> \param grid_atom ...
220 : !> \param ik ...
221 : ! **************************************************************************************************
222 4280 : SUBROUTINE calculate_g0(rho0_mpole, grid_atom, ik)
223 :
224 : TYPE(rho0_mpole_type), POINTER :: rho0_mpole
225 : TYPE(grid_atom_type), POINTER :: grid_atom
226 : INTEGER, INTENT(IN) :: ik
227 :
228 : INTEGER :: ir, l, lmax, nr
229 : REAL(dp) :: c1, prefactor, root_z_h, z_h
230 4280 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: erf_z_h, gexp, gh_tmp, int1, int2
231 :
232 4280 : nr = grid_atom%nr
233 4280 : lmax = rho0_mpole%lmax0_kind(ik)
234 4280 : z_h = rho0_mpole%zet0_h
235 4280 : root_z_h = SQRT(z_h)
236 :
237 : ! Allocate g0
238 4280 : CALL reallocate(rho0_mpole%mp_gau(ik)%g0_h, 1, nr, 0, lmax)
239 4280 : CALL reallocate(rho0_mpole%mp_gau(ik)%vg0_h, 1, nr, 0, lmax)
240 :
241 29960 : ALLOCATE (gexp(nr), gh_tmp(nr), erf_z_h(nr), int1(nr), int2(nr))
242 :
243 223200 : gh_tmp(1:nr) = EXP(-z_h*grid_atom%rad2(1:nr))
244 :
245 223200 : DO ir = 1, nr
246 223200 : erf_z_h(ir) = erf(grid_atom%rad(ir)*root_z_h)
247 : END DO
248 :
249 223200 : DO ir = 1, nr
250 223200 : IF (gh_tmp(ir) < 1.0E-30_dp) gh_tmp(ir) = 0.0_dp
251 : END DO
252 :
253 223200 : gexp(1:nr) = gh_tmp(1:nr)
254 : rho0_mpole%mp_gau(ik)%g0_h(1:nr, 0) = gh_tmp(1:nr)* &
255 223200 : rho0_mpole%norm_g0l_h(0)
256 4280 : CALL whittaker_c0a(int1, grid_atom%rad, gh_tmp, erf_z_h, z_h, 0, 0, nr)
257 4280 : CALL whittaker_ci(int2, grid_atom%rad, gh_tmp, z_h, 0, nr)
258 :
259 4280 : prefactor = fourpi*rho0_mpole%norm_g0l_h(0)
260 :
261 4280 : c1 = SQRT(pi*pi*pi/(z_h*z_h*z_h))*rho0_mpole%norm_g0l_h(0)
262 :
263 223200 : DO ir = 1, nr
264 223200 : rho0_mpole%mp_gau(ik)%vg0_h(ir, 0) = c1*erf_z_h(ir)*grid_atom%oorad2l(ir, 1)
265 : END DO
266 :
267 11724 : DO l = 1, lmax
268 388604 : gh_tmp(1:nr) = gh_tmp(1:nr)*grid_atom%rad(1:nr)
269 : rho0_mpole%mp_gau(ik)%g0_h(1:nr, l) = gh_tmp(1:nr)* &
270 388604 : rho0_mpole%norm_g0l_h(l)
271 :
272 7444 : prefactor = fourpi/(2.0_dp*l + 1.0_dp)*rho0_mpole%norm_g0l_h(l)
273 7444 : CALL whittaker_c0a(int1, grid_atom%rad, gexp, erf_z_h, z_h, l, l, nr)
274 392884 : DO ir = 1, nr
275 : rho0_mpole%mp_gau(ik)%vg0_h(ir, l) = prefactor*(int1(ir) + &
276 388604 : int2(ir)*grid_atom%rad2l(ir, l))
277 : END DO
278 :
279 : END DO ! l
280 :
281 4280 : DEALLOCATE (gexp, erf_z_h, gh_tmp, int1, int2)
282 4280 : END SUBROUTINE calculate_g0
283 :
284 : ! **************************************************************************************************
285 : !> \brief ...
286 : !> \param mp_gau ...
287 : ! **************************************************************************************************
288 2168 : SUBROUTINE deallocate_mpole_gau(mp_gau)
289 :
290 : TYPE(mpole_gau_overlap), DIMENSION(:), POINTER :: mp_gau
291 :
292 : INTEGER :: ikind, nkind
293 :
294 2168 : nkind = SIZE(mp_gau)
295 :
296 6448 : DO ikind = 1, nkind
297 :
298 4280 : IF (ASSOCIATED(mp_gau(ikind)%Qlm_gg)) THEN
299 3926 : DEALLOCATE (mp_gau(ikind)%Qlm_gg)
300 : END IF
301 :
302 4280 : DEALLOCATE (mp_gau(ikind)%g0_h)
303 :
304 6448 : DEALLOCATE (mp_gau(ikind)%vg0_h)
305 : END DO
306 :
307 2168 : DEALLOCATE (mp_gau)
308 :
309 2168 : END SUBROUTINE deallocate_mpole_gau
310 :
311 : ! **************************************************************************************************
312 : !> \brief ...
313 : !> \param mp_rho ...
314 : ! **************************************************************************************************
315 2168 : SUBROUTINE deallocate_mpole_rho(mp_rho)
316 :
317 : TYPE(mpole_rho_atom), DIMENSION(:), POINTER :: mp_rho
318 :
319 : INTEGER :: iat, natom
320 :
321 2168 : natom = SIZE(mp_rho)
322 :
323 9216 : DO iat = 1, natom
324 7048 : DEALLOCATE (mp_rho(iat)%Qlm_h)
325 7048 : DEALLOCATE (mp_rho(iat)%Qlm_s)
326 7048 : DEALLOCATE (mp_rho(iat)%Qlm_tot)
327 9216 : DEALLOCATE (mp_rho(iat)%Qlm_car)
328 : END DO
329 :
330 2168 : DEALLOCATE (mp_rho)
331 :
332 2168 : END SUBROUTINE deallocate_mpole_rho
333 :
334 : ! **************************************************************************************************
335 : !> \brief ...
336 : !> \param rho0_atom_set ...
337 : ! **************************************************************************************************
338 2168 : SUBROUTINE deallocate_rho0_atom(rho0_atom_set)
339 :
340 : TYPE(rho0_atom_type), DIMENSION(:), POINTER :: rho0_atom_set
341 :
342 : INTEGER :: iat, natom
343 :
344 2168 : IF (ASSOCIATED(rho0_atom_set)) THEN
345 :
346 2168 : natom = SIZE(rho0_atom_set)
347 :
348 9216 : DO iat = 1, natom
349 7048 : IF (ASSOCIATED(rho0_atom_set(iat)%rho0_rad_h)) THEN
350 7048 : DEALLOCATE (rho0_atom_set(iat)%rho0_rad_h%r_coef)
351 7048 : DEALLOCATE (rho0_atom_set(iat)%rho0_rad_h)
352 : END IF
353 9216 : IF (ASSOCIATED(rho0_atom_set(iat)%vrho0_rad_h)) THEN
354 7048 : DEALLOCATE (rho0_atom_set(iat)%vrho0_rad_h%r_coef)
355 7048 : DEALLOCATE (rho0_atom_set(iat)%vrho0_rad_h)
356 : END IF
357 : END DO
358 :
359 2168 : DEALLOCATE (rho0_atom_set)
360 : ELSE
361 : CALL cp_abort(__LOCATION__, &
362 : "The pointer rho0_atom_set is not associated and "// &
363 0 : "cannot be deallocated")
364 : END IF
365 :
366 2168 : END SUBROUTINE deallocate_rho0_atom
367 : ! **************************************************************************************************
368 : !> \brief ...
369 : !> \param rho0 ...
370 : ! **************************************************************************************************
371 2168 : SUBROUTINE deallocate_rho0_mpole(rho0)
372 :
373 : TYPE(rho0_mpole_type), POINTER :: rho0
374 :
375 2168 : IF (ASSOCIATED(rho0)) THEN
376 :
377 2168 : IF (ASSOCIATED(rho0%mp_gau)) CALL deallocate_mpole_gau(rho0%mp_gau)
378 :
379 2168 : IF (ASSOCIATED(rho0%mp_rho)) CALL deallocate_mpole_rho(rho0%mp_rho)
380 :
381 2168 : IF (ASSOCIATED(rho0%lmax0_kind)) THEN
382 2168 : DEALLOCATE (rho0%lmax0_kind)
383 : END IF
384 :
385 2168 : IF (ASSOCIATED(rho0%norm_g0l_h)) THEN
386 2168 : DEALLOCATE (rho0%norm_g0l_h)
387 : END IF
388 :
389 2168 : IF (ASSOCIATED(rho0%rho0_s_rs)) THEN
390 2168 : CALL rho0%rho0_s_rs%release()
391 2168 : DEALLOCATE (rho0%rho0_s_rs)
392 : END IF
393 :
394 2168 : IF (ASSOCIATED(rho0%rho0_s_gs)) THEN
395 2168 : CALL rho0%rho0_s_gs%release()
396 2168 : DEALLOCATE (rho0%rho0_s_gs)
397 : END IF
398 :
399 2168 : IF (ASSOCIATED(rho0%rhoz_cneo_s_rs)) THEN
400 8 : CALL rho0%rhoz_cneo_s_rs%release()
401 8 : DEALLOCATE (rho0%rhoz_cneo_s_rs)
402 : END IF
403 :
404 2168 : IF (ASSOCIATED(rho0%rhoz_cneo_s_gs)) THEN
405 8 : CALL rho0%rhoz_cneo_s_gs%release()
406 8 : DEALLOCATE (rho0%rhoz_cneo_s_gs)
407 : END IF
408 :
409 2168 : DEALLOCATE (rho0)
410 : ELSE
411 : CALL cp_abort(__LOCATION__, &
412 : "The pointer rho0 is not associated and "// &
413 0 : "cannot be deallocated")
414 : END IF
415 :
416 2168 : END SUBROUTINE deallocate_rho0_mpole
417 :
418 : ! **************************************************************************************************
419 : !> \brief ...
420 : !> \param rho0_mpole ...
421 : !> \param g0_h ...
422 : !> \param vg0_h ...
423 : !> \param iat ...
424 : !> \param ikind ...
425 : !> \param lmax_0 ...
426 : !> \param l0_ikind ...
427 : !> \param mp_gau_ikind ...
428 : !> \param mp_rho ...
429 : !> \param norm_g0l_h ...
430 : !> \param Qlm_gg ...
431 : !> \param Qlm_car ...
432 : !> \param Qlm_tot ...
433 : !> \param zet0_h ...
434 : !> \param igrid_zet0_s ...
435 : !> \param rpgf0_h ...
436 : !> \param rpgf0_s ...
437 : !> \param max_rpgf0_s ...
438 : !> \param rho0_s_rs ...
439 : !> \param rho0_s_gs ...
440 : !> \param rhoz_cneo_s_rs ...
441 : !> \param rhoz_cneo_s_gs ...
442 : ! **************************************************************************************************
443 247486 : SUBROUTINE get_rho0_mpole(rho0_mpole, g0_h, vg0_h, iat, ikind, lmax_0, l0_ikind, &
444 : mp_gau_ikind, mp_rho, norm_g0l_h, &
445 : Qlm_gg, Qlm_car, Qlm_tot, &
446 : zet0_h, igrid_zet0_s, rpgf0_h, rpgf0_s, &
447 : max_rpgf0_s, rho0_s_rs, rho0_s_gs, &
448 : rhoz_cneo_s_rs, rhoz_cneo_s_gs)
449 :
450 : TYPE(rho0_mpole_type), POINTER :: rho0_mpole
451 : REAL(dp), DIMENSION(:, :), OPTIONAL, POINTER :: g0_h, vg0_h
452 : INTEGER, INTENT(IN), OPTIONAL :: iat, ikind
453 : INTEGER, INTENT(OUT), OPTIONAL :: lmax_0, l0_ikind
454 : TYPE(mpole_gau_overlap), OPTIONAL, POINTER :: mp_gau_ikind
455 : TYPE(mpole_rho_atom), DIMENSION(:), OPTIONAL, &
456 : POINTER :: mp_rho
457 : REAL(dp), DIMENSION(:), OPTIONAL, POINTER :: norm_g0l_h
458 : REAL(dp), DIMENSION(:, :, :), OPTIONAL, POINTER :: Qlm_gg
459 : REAL(dp), DIMENSION(:), OPTIONAL, POINTER :: Qlm_car, Qlm_tot
460 : REAL(dp), INTENT(OUT), OPTIONAL :: zet0_h
461 : INTEGER, INTENT(OUT), OPTIONAL :: igrid_zet0_s
462 : REAL(dp), INTENT(OUT), OPTIONAL :: rpgf0_h, rpgf0_s, max_rpgf0_s
463 : TYPE(pw_r3d_rs_type), OPTIONAL, POINTER :: rho0_s_rs
464 : TYPE(pw_c1d_gs_type), OPTIONAL, POINTER :: rho0_s_gs
465 : TYPE(pw_r3d_rs_type), OPTIONAL, POINTER :: rhoz_cneo_s_rs
466 : TYPE(pw_c1d_gs_type), OPTIONAL, POINTER :: rhoz_cneo_s_gs
467 :
468 247486 : IF (ASSOCIATED(rho0_mpole)) THEN
469 :
470 247486 : IF (PRESENT(lmax_0)) lmax_0 = rho0_mpole%lmax_0
471 247486 : IF (PRESENT(mp_rho)) mp_rho => rho0_mpole%mp_rho
472 247486 : IF (PRESENT(norm_g0l_h)) norm_g0l_h => rho0_mpole%norm_g0l_h
473 247486 : IF (PRESENT(zet0_h)) zet0_h = rho0_mpole%zet0_h
474 247486 : IF (PRESENT(igrid_zet0_s)) igrid_zet0_s = rho0_mpole%igrid_zet0_s
475 247486 : IF (PRESENT(max_rpgf0_s)) max_rpgf0_s = rho0_mpole%max_rpgf0_s
476 247486 : IF (PRESENT(rho0_s_rs)) rho0_s_rs => rho0_mpole%rho0_s_rs
477 247486 : IF (PRESENT(rho0_s_gs)) rho0_s_gs => rho0_mpole%rho0_s_gs
478 247486 : IF (PRESENT(rhoz_cneo_s_rs)) rhoz_cneo_s_rs => rho0_mpole%rhoz_cneo_s_rs
479 247486 : IF (PRESENT(rhoz_cneo_s_gs)) rhoz_cneo_s_gs => rho0_mpole%rhoz_cneo_s_gs
480 :
481 247486 : IF (PRESENT(ikind)) THEN
482 139280 : IF (PRESENT(l0_ikind)) l0_ikind = rho0_mpole%lmax0_kind(ikind)
483 139280 : IF (PRESENT(mp_gau_ikind)) mp_gau_ikind => rho0_mpole%mp_gau(ikind)
484 139280 : IF (PRESENT(g0_h)) g0_h => rho0_mpole%mp_gau(ikind)%g0_h
485 139280 : IF (PRESENT(vg0_h)) vg0_h => rho0_mpole%mp_gau(ikind)%vg0_h
486 139280 : IF (PRESENT(Qlm_gg)) Qlm_gg => rho0_mpole%mp_gau(ikind)%Qlm_gg
487 139280 : IF (PRESENT(rpgf0_h)) rpgf0_h = rho0_mpole%mp_gau(ikind)%rpgf0_h
488 139280 : IF (PRESENT(rpgf0_s)) rpgf0_s = rho0_mpole%mp_gau(ikind)%rpgf0_s
489 : END IF
490 247486 : IF (PRESENT(iat)) THEN
491 53820 : IF (PRESENT(Qlm_car)) Qlm_car => rho0_mpole%mp_rho(iat)%Qlm_car
492 53820 : IF (PRESENT(Qlm_tot)) Qlm_tot => rho0_mpole%mp_rho(iat)%Qlm_tot
493 : END IF
494 :
495 : ELSE
496 0 : CPABORT("The pointer rho0_mpole is not associated")
497 : END IF
498 :
499 247486 : END SUBROUTINE get_rho0_mpole
500 :
501 : ! **************************************************************************************************
502 : !> \brief ...
503 : !> \param mp_rho ...
504 : !> \param nchan_s ...
505 : !> \param nchan_c ...
506 : !> \param zeff ...
507 : ! **************************************************************************************************
508 7048 : SUBROUTINE initialize_mpole_rho(mp_rho, nchan_s, nchan_c, zeff)
509 :
510 : TYPE(mpole_rho_atom) :: mp_rho
511 : INTEGER, INTENT(IN) :: nchan_s, nchan_c
512 : REAL(KIND=dp), INTENT(IN) :: zeff
513 :
514 7048 : CALL reallocate(mp_rho%Qlm_h, 1, nchan_s)
515 7048 : CALL reallocate(mp_rho%Qlm_s, 1, nchan_s)
516 7048 : CALL reallocate(mp_rho%Qlm_tot, 1, nchan_s)
517 7048 : CALL reallocate(mp_rho%Qlm_car, 1, nchan_c)
518 :
519 60176 : mp_rho%Qlm_h = 0.0_dp
520 60176 : mp_rho%Qlm_s = 0.0_dp
521 60176 : mp_rho%Qlm_tot = 0.0_dp
522 66804 : mp_rho%Qlm_car = 0.0_dp
523 7048 : mp_rho%Qlm_z = -2.0_dp*rootpi*Zeff
524 21144 : mp_rho%Q0 = 0.0_dp
525 :
526 7048 : END SUBROUTINE initialize_mpole_rho
527 :
528 : ! **************************************************************************************************
529 : !> \brief ...
530 : !> \param rho0_mpole ...
531 : !> \param unit_str ...
532 : !> \param output_unit ...
533 : ! **************************************************************************************************
534 1 : SUBROUTINE write_rho0_info(rho0_mpole, unit_str, output_unit)
535 :
536 : TYPE(rho0_mpole_type), POINTER :: rho0_mpole
537 : CHARACTER(LEN=*), INTENT(IN) :: unit_str
538 : INTEGER, INTENT(in) :: output_unit
539 :
540 : INTEGER :: ikind, l, nkind
541 : REAL(dp) :: conv
542 :
543 1 : IF (ASSOCIATED(rho0_mpole)) THEN
544 1 : conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
545 :
546 : WRITE (UNIT=output_unit, FMT="(/,T5,A,/)") &
547 1 : "*** Compensation density charges data set ***"
548 : WRITE (UNIT=output_unit, FMT="(T2,A,T35,f16.10)") &
549 1 : "- Rho0 exponent :", rho0_mpole%zet0_h
550 : WRITE (UNIT=output_unit, FMT="(T2,A,T35,I5)") &
551 1 : "- Global max l :", rho0_mpole%lmax_0
552 :
553 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
554 1 : "- Normalization constants for g0"
555 2 : DO l = 0, rho0_mpole%lmax_0
556 : WRITE (UNIT=output_unit, FMT="(T20,A,T31,I2,T38,A,f15.5)") &
557 2 : "ang. mom.= ", l, " hard= ", rho0_mpole%norm_g0l_h(l)
558 : END DO
559 :
560 1 : nkind = SIZE(rho0_mpole%lmax0_kind, 1)
561 2 : DO ikind = 1, nkind
562 : WRITE (UNIT=output_unit, FMT="(/,T2,A,T55,I2)") &
563 : "- rho0 max L and radii in "//TRIM(unit_str)// &
564 1 : " for the atom kind :", ikind
565 :
566 : WRITE (UNIT=output_unit, FMT="(T2,T20,A,T55,I5)") &
567 1 : "=> l max :", rho0_mpole%lmax0_kind(ikind)
568 :
569 : WRITE (UNIT=output_unit, FMT="(T2,T20,A,T55,f20.10)") &
570 1 : "=> max radius of g0: ", &
571 3 : rho0_mpole%mp_gau(ikind)%rpgf0_h*conv
572 : END DO ! ikind
573 :
574 : ELSE
575 : WRITE (UNIT=output_unit, FMT="(/,T5,A,/)") &
576 0 : ' WARNING: I cannot print rho0, it is not associated'
577 : END IF
578 :
579 1 : END SUBROUTINE write_rho0_info
580 0 : END MODULE qs_rho0_types
|
