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 : ! **************************************************************************************************
9 : MODULE qs_rho0_methods
10 :
11 : USE ao_util, ONLY: exp_radius,&
12 : gaussint_sph,&
13 : trace_r_AxB
14 : USE atomic_kind_types, ONLY: atomic_kind_type,&
15 : get_atomic_kind,&
16 : get_atomic_kind_set
17 : USE basis_set_types, ONLY: get_gto_basis_set,&
18 : gto_basis_set_type
19 : USE cp_control_types, ONLY: gapw_control_type
20 : USE cp_log_handling, ONLY: cp_get_default_logger,&
21 : cp_logger_type
22 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
23 : cp_print_key_unit_nr
24 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
25 : section_vals_type,&
26 : section_vals_val_get
27 : USE kinds, ONLY: default_string_length,&
28 : dp
29 : USE mathconstants, ONLY: fourpi
30 : USE memory_utilities, ONLY: reallocate
31 : USE orbital_pointers, ONLY: indco,&
32 : indso,&
33 : nco,&
34 : ncoset,&
35 : nso,&
36 : nsoset
37 : USE orbital_transformation_matrices, ONLY: orbtramat
38 : USE qs_cneo_methods, ONLY: allocate_rhoz_cneo_internals,&
39 : init_cneo_potential_internals
40 : USE qs_cneo_types, ONLY: cneo_potential_type,&
41 : rhoz_cneo_type
42 : USE qs_cneo_utils, ONLY: cneo_scatter
43 : USE qs_environment_types, ONLY: get_qs_env,&
44 : qs_environment_type
45 : USE qs_grid_atom, ONLY: grid_atom_type
46 : USE qs_harmonics_atom, ONLY: get_none0_cg_list,&
47 : harmonics_atom_type
48 : USE qs_kind_types, ONLY: get_qs_kind,&
49 : get_qs_kind_set,&
50 : qs_kind_type,&
51 : set_qs_kind
52 : USE qs_local_rho_types, ONLY: allocate_rhoz,&
53 : calculate_rhoz,&
54 : local_rho_type,&
55 : rhoz_type,&
56 : set_local_rho
57 : USE qs_oce_methods, ONLY: prj_scatter
58 : USE qs_rho0_types, ONLY: &
59 : allocate_multipoles, allocate_rho0_atom, allocate_rho0_atom_rad, allocate_rho0_mpole, &
60 : calculate_g0, get_rho0_mpole, initialize_mpole_rho, mpole_gau_overlap, mpole_rho_atom, &
61 : rho0_atom_type, rho0_mpole_type, write_rho0_info
62 : USE qs_rho_atom_types, ONLY: get_rho_atom,&
63 : rho_atom_coeff,&
64 : rho_atom_type
65 : #include "./base/base_uses.f90"
66 :
67 : IMPLICIT NONE
68 :
69 : PRIVATE
70 :
71 : ! Global parameters (only in this module)
72 :
73 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_rho0_methods'
74 :
75 : ! Public subroutines
76 :
77 : PUBLIC :: calculate_rho0_atom, init_rho0
78 :
79 : CONTAINS
80 :
81 : ! **************************************************************************************************
82 : !> \brief ...
83 : !> \param Qlm_gg ...
84 : !> \param basis_1c ...
85 : !> \param harmonics ...
86 : !> \param nchannels ...
87 : !> \param nsotot ...
88 : ! **************************************************************************************************
89 3934 : SUBROUTINE calculate_mpole_gau(Qlm_gg, basis_1c, harmonics, nchannels, nsotot)
90 :
91 : REAL(dp), DIMENSION(:, :, :), POINTER :: Qlm_gg
92 : TYPE(gto_basis_set_type), POINTER :: basis_1c
93 : TYPE(harmonics_atom_type), POINTER :: harmonics
94 : INTEGER, INTENT(IN) :: nchannels, nsotot
95 :
96 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_mpole_gau'
97 :
98 : INTEGER :: handle, icg, ig1, ig2, ipgf1, ipgf2, iset1, iset2, iso, iso1, iso2, l, l1, l2, &
99 : llmax, m1, m2, max_iso_not0_local, max_s_harm, maxl, maxso, n1, n2, nset
100 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list
101 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list
102 3934 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf
103 : REAL(KIND=dp) :: zet1, zet2
104 3934 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zet
105 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: my_CG
106 :
107 3934 : CALL timeset(routineN, handle)
108 :
109 3934 : NULLIFY (lmax, lmin, npgf, my_CG, zet)
110 :
111 : CALL reallocate(Qlm_gg, 1, nsotot, 1, nsotot, 1, &
112 3934 : MIN(nchannels, harmonics%max_iso_not0))
113 :
114 : CALL get_gto_basis_set(gto_basis_set=basis_1c, &
115 : lmax=lmax, lmin=lmin, maxso=maxso, &
116 3934 : npgf=npgf, nset=nset, zet=zet, maxl=maxl)
117 :
118 3934 : max_s_harm = harmonics%max_s_harm
119 3934 : llmax = harmonics%llmax
120 :
121 23604 : ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm))
122 :
123 3934 : my_CG => harmonics%my_CG
124 :
125 3934 : m1 = 0
126 15084 : DO iset1 = 1, nset
127 : m2 = 0
128 52888 : DO iset2 = 1, nset
129 :
130 : CALL get_none0_cg_list(my_CG, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
131 41738 : max_s_harm, llmax, cg_list, cg_n_list, max_iso_not0_local)
132 :
133 41738 : n1 = nsoset(lmax(iset1))
134 138514 : DO ipgf1 = 1, npgf(iset1)
135 96776 : zet1 = zet(ipgf1, iset1)
136 :
137 96776 : n2 = nsoset(lmax(iset2))
138 391372 : DO ipgf2 = 1, npgf(iset2)
139 252858 : zet2 = zet(ipgf2, iset2)
140 :
141 1799780 : DO iso = 1, MIN(nchannels, max_iso_not0_local)
142 1450146 : l = indso(1, iso)
143 4286378 : DO icg = 1, cg_n_list(iso)
144 2583374 : iso1 = cg_list(1, icg, iso)
145 2583374 : iso2 = cg_list(2, icg, iso)
146 :
147 2583374 : l1 = indso(1, iso1)
148 2583374 : l2 = indso(1, iso2)
149 2583374 : ig1 = iso1 + n1*(ipgf1 - 1) + m1
150 2583374 : ig2 = iso2 + n2*(ipgf2 - 1) + m2
151 :
152 : Qlm_gg(ig1, ig2, iso) = fourpi/(2._dp*l + 1._dp)* &
153 4033520 : my_CG(iso1, iso2, iso)*gaussint_sph(zet1 + zet2, l + l1 + l2)
154 : END DO ! icg
155 : END DO ! iso
156 :
157 : END DO ! ipgf2
158 : END DO ! ipgf1
159 94626 : m2 = m2 + maxso
160 : END DO ! iset2
161 15084 : m1 = m1 + maxso
162 : END DO ! iset1
163 :
164 3934 : DEALLOCATE (cg_list, cg_n_list)
165 :
166 3934 : CALL timestop(handle)
167 7868 : END SUBROUTINE calculate_mpole_gau
168 :
169 : ! **************************************************************************************************
170 : !> \brief ...
171 : !> \param gapw_control ...
172 : !> \param rho_atom_set ...
173 : !> \param rhoz_cneo_set ...
174 : !> \param rho0_atom_set ...
175 : !> \param rho0_mp ...
176 : !> \param a_list ...
177 : !> \param natom ...
178 : !> \param ikind ...
179 : !> \param qs_kind ...
180 : !> \param rho0_h_tot ...
181 : ! **************************************************************************************************
182 36308 : SUBROUTINE calculate_rho0_atom(gapw_control, rho_atom_set, rhoz_cneo_set, rho0_atom_set, &
183 36308 : rho0_mp, a_list, natom, ikind, qs_kind, rho0_h_tot)
184 :
185 : TYPE(gapw_control_type), POINTER :: gapw_control
186 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
187 : TYPE(rhoz_cneo_type), DIMENSION(:), POINTER :: rhoz_cneo_set
188 : TYPE(rho0_atom_type), DIMENSION(:), POINTER :: rho0_atom_set
189 : TYPE(rho0_mpole_type), POINTER :: rho0_mp
190 : INTEGER, DIMENSION(:), INTENT(IN) :: a_list
191 : INTEGER, INTENT(IN) :: natom, ikind
192 : TYPE(qs_kind_type), INTENT(IN) :: qs_kind
193 : REAL(KIND=dp), INTENT(INOUT) :: rho0_h_tot
194 :
195 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_rho0_atom'
196 :
197 : INTEGER :: handle, iat, iatom, ic, ico, ir, is, &
198 : iso, ispin, l, lmax0, lshell, lx, ly, &
199 : lz, nr, nsotot, nsotot_nuc, nspins
200 : LOGICAL :: paw_atom
201 : REAL(KIND=dp) :: sum1
202 36308 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cpc_h_nuc, cpc_s_nuc
203 36308 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: cpc_ah, cpc_as
204 36308 : REAL(KIND=dp), DIMENSION(:), POINTER :: norm_g0l_h
205 36308 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: g0_h, vg0_h
206 : TYPE(cneo_potential_type), POINTER :: cneo_potential
207 : TYPE(grid_atom_type), POINTER :: g_atom
208 : TYPE(harmonics_atom_type), POINTER :: harmonics
209 : TYPE(mpole_gau_overlap), POINTER :: mpole_gau
210 : TYPE(mpole_rho_atom), POINTER :: mpole_rho
211 36308 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: cpc_h, cpc_s
212 : TYPE(rho_atom_type), POINTER :: rho_atom
213 :
214 36308 : CALL timeset(routineN, handle)
215 :
216 36308 : NULLIFY (mpole_gau)
217 36308 : NULLIFY (mpole_rho)
218 36308 : NULLIFY (g0_h, vg0_h, g_atom)
219 36308 : NULLIFY (norm_g0l_h, harmonics)
220 36308 : NULLIFY (cneo_potential)
221 :
222 : CALL get_rho0_mpole(rho0_mpole=rho0_mp, ikind=ikind, &
223 : l0_ikind=lmax0, mp_gau_ikind=mpole_gau, &
224 : g0_h=g0_h, &
225 : vg0_h=vg0_h, &
226 36308 : norm_g0l_h=norm_g0l_h)
227 :
228 : CALL get_qs_kind(qs_kind, harmonics=harmonics, paw_atom=paw_atom, grid_atom=g_atom, &
229 36308 : cneo_potential=cneo_potential)
230 :
231 36308 : nr = g_atom%nr
232 :
233 : ! Set density coefficient to zero before the calculation
234 96930 : DO iat = 1, natom
235 60622 : iatom = a_list(iat)
236 23233856 : rho0_atom_set(iatom)%rho0_rad_h%r_coef = 0.0_dp
237 483276 : rho0_mp%mp_rho(iatom)%Qlm_tot = 0.0_dp
238 : ! When no nuclear density matrix is available, use spherical zeff
239 60622 : IF (.NOT. ASSOCIATED(cneo_potential)) THEN
240 60530 : rho0_mp%mp_rho(iatom)%Qlm_tot(1) = rho0_mp%mp_rho(iatom)%Qlm_z
241 : ELSE
242 92 : IF (.NOT. rhoz_cneo_set(iatom)%ready) THEN
243 14 : rho0_mp%mp_rho(iatom)%Qlm_tot(1) = rho0_mp%mp_rho(iatom)%Qlm_z
244 : END IF
245 : END IF
246 181866 : rho0_mp%mp_rho(iatom)%Q0 = 0.0_dp
247 568646 : rho0_mp%mp_rho(iatom)%Qlm_car = 0.0_dp
248 : END DO
249 :
250 36308 : IF (.NOT. (.NOT. paw_atom .AND. gapw_control%nopaw_as_gpw)) THEN
251 90532 : DO iat = 1, natom
252 55864 : iatom = a_list(iat)
253 55864 : mpole_rho => rho0_mp%mp_rho(iatom)
254 55864 : rho_atom => rho_atom_set(iatom)
255 :
256 55864 : IF (paw_atom) THEN
257 55864 : NULLIFY (cpc_h, cpc_s)
258 55864 : CALL get_rho_atom(rho_atom=rho_atom, cpc_h=cpc_h, cpc_s=cpc_s)
259 55864 : nspins = SIZE(cpc_h)
260 55864 : nsotot = SIZE(mpole_gau%Qlm_gg, 1)
261 279320 : ALLOCATE (cpc_ah(nsotot, nsotot, nspins))
262 175530360 : cpc_ah = 0._dp
263 223456 : ALLOCATE (cpc_as(nsotot, nsotot, nspins))
264 175530360 : cpc_as = 0._dp
265 121016 : DO ispin = 1, nspins
266 65152 : CALL prj_scatter(cpc_h(ispin)%r_coef, cpc_ah(:, :, ispin), qs_kind)
267 121016 : CALL prj_scatter(cpc_s(ispin)%r_coef, cpc_as(:, :, ispin), qs_kind)
268 : END DO
269 55864 : nsotot_nuc = 0
270 55864 : IF (ASSOCIATED(cneo_potential)) THEN
271 92 : IF (rhoz_cneo_set(iatom)%ready) THEN
272 78 : nsotot_nuc = SIZE(cneo_potential%Qlm_gg, 1)
273 468 : ALLOCATE (cpc_h_nuc(nsotot_nuc, nsotot_nuc), cpc_s_nuc(nsotot_nuc, nsotot_nuc))
274 518154 : cpc_h_nuc = 0._dp
275 518154 : cpc_s_nuc = 0._dp
276 : CALL cneo_scatter(rhoz_cneo_set(iatom)%cpc_h, cpc_h_nuc, cneo_potential%npsgf, &
277 78 : cneo_potential%n2oindex)
278 : CALL cneo_scatter(rhoz_cneo_set(iatom)%cpc_s, cpc_s_nuc, cneo_potential%npsgf, &
279 78 : cneo_potential%n2oindex)
280 : END IF
281 : END IF
282 :
283 : ! Total charge (hard-soft) at atom
284 : IF (paw_atom) THEN
285 121016 : DO ispin = 1, nspins
286 : mpole_rho%Q0(ispin) = (trace_r_AxB(mpole_gau%Qlm_gg(:, :, 1), nsotot, &
287 : cpc_ah(:, :, ispin), nsotot, nsotot, nsotot) &
288 : - trace_r_AxB(mpole_gau%Qlm_gg(:, :, 1), nsotot, &
289 121016 : cpc_as(:, :, ispin), nsotot, nsotot, nsotot))/SQRT(fourpi)
290 : END DO
291 : END IF
292 : ! Multipoles of local charge distribution
293 473760 : DO iso = 1, MIN(nsoset(lmax0), harmonics%max_iso_not0)
294 55864 : IF (paw_atom) THEN
295 417896 : mpole_rho%Qlm_h(iso) = 0.0_dp
296 417896 : mpole_rho%Qlm_s(iso) = 0.0_dp
297 :
298 900408 : DO ispin = 1, nspins
299 : mpole_rho%Qlm_h(iso) = mpole_rho%Qlm_h(iso) + &
300 : trace_r_AxB(mpole_gau%Qlm_gg(:, :, iso), nsotot, &
301 482512 : cpc_ah(:, :, ispin), nsotot, nsotot, nsotot)
302 : mpole_rho%Qlm_s(iso) = mpole_rho%Qlm_s(iso) + &
303 : trace_r_AxB(mpole_gau%Qlm_gg(:, :, iso), nsotot, &
304 900408 : cpc_as(:, :, ispin), nsotot, nsotot, nsotot)
305 : END DO ! ispin
306 :
307 : mpole_rho%Qlm_tot(iso) = mpole_rho%Qlm_tot(iso) + &
308 417896 : mpole_rho%Qlm_h(iso) - mpole_rho%Qlm_s(iso)
309 : END IF
310 : END DO ! iso
311 :
312 : ! Multipoles of CNEO quantum nuclear charge distribuition
313 55864 : IF (ASSOCIATED(cneo_potential)) THEN
314 92 : IF (rhoz_cneo_set(iatom)%ready) THEN
315 780 : DO iso = 1, MIN(nsoset(lmax0), cneo_potential%harmonics%max_iso_not0)
316 : mpole_rho%Qlm_tot(iso) = mpole_rho%Qlm_tot(iso) - cneo_potential%zeff* &
317 : trace_r_AxB(cneo_potential%Qlm_gg(:, :, iso), nsotot_nuc, &
318 : cpc_h_nuc - cpc_s_nuc, nsotot_nuc, &
319 4663464 : nsotot_nuc, nsotot_nuc)
320 : END DO ! iso
321 : END IF
322 : END IF
323 :
324 55864 : DEALLOCATE (cpc_ah, cpc_as)
325 55864 : IF (ALLOCATED(cpc_h_nuc)) DEALLOCATE (cpc_h_nuc)
326 55864 : IF (ALLOCATED(cpc_s_nuc)) DEALLOCATE (cpc_s_nuc)
327 : END IF
328 :
329 510068 : DO iso = 1, nsoset(lmax0)
330 417896 : l = indso(1, iso)
331 : rho0_atom_set(iatom)%rho0_rad_h%r_coef(1:nr, iso) = &
332 45443736 : g0_h(1:nr, l)*mpole_rho%Qlm_tot(iso)
333 : rho0_atom_set(iatom)%vrho0_rad_h%r_coef(1:nr, iso) = &
334 45443736 : vg0_h(1:nr, l)*mpole_rho%Qlm_tot(iso)
335 :
336 : ! When CNEO is enabled, it is possible for rho0 to have a higher angualr momentum
337 : ! than that of the electronic density. In that case, the nuclear density must have
338 : ! a higher angular momentum, but cneo_potential%harmonics%slm_int is not initialized.
339 : ! For higher angular momenta, simply use the fact that slm_int(iso>1)=0
340 473760 : IF (iso <= harmonics%max_iso_not0) THEN
341 : sum1 = 0.0_dp
342 22930816 : DO ir = 1, nr
343 : sum1 = sum1 + g_atom%wr(ir)* &
344 22930816 : rho0_atom_set(iatom)%rho0_rad_h%r_coef(ir, iso)
345 : END DO
346 417896 : rho0_h_tot = rho0_h_tot + sum1*harmonics%slm_int(iso)
347 : END IF
348 : END DO ! iso
349 : END DO ! iat
350 : END IF
351 :
352 : ! Transform the coefficinets from spherical to Cartesian
353 36308 : IF (.NOT. paw_atom .AND. gapw_control%nopaw_as_gpw) THEN
354 6398 : DO iat = 1, natom
355 4758 : iatom = a_list(iat)
356 4758 : mpole_rho => rho0_mp%mp_rho(iatom)
357 :
358 11156 : DO lshell = 0, lmax0
359 14274 : DO ic = 1, nco(lshell)
360 4758 : ico = ic + ncoset(lshell - 1)
361 9516 : mpole_rho%Qlm_car(ico) = 0.0_dp
362 : END DO
363 : END DO
364 : END DO
365 : ELSE
366 90532 : DO iat = 1, natom
367 55864 : iatom = a_list(iat)
368 55864 : mpole_rho => rho0_mp%mp_rho(iatom)
369 235816 : DO lshell = 0, lmax0
370 668106 : DO ic = 1, nco(lshell)
371 466958 : ico = ic + ncoset(lshell - 1)
372 466958 : mpole_rho%Qlm_car(ico) = 0.0_dp
373 466958 : lx = indco(1, ico)
374 466958 : ly = indco(2, ico)
375 466958 : lz = indco(3, ico)
376 2502100 : DO is = 1, nso(lshell)
377 1889858 : iso = is + nsoset(lshell - 1)
378 : mpole_rho%Qlm_car(ico) = mpole_rho%Qlm_car(ico) + &
379 : norm_g0l_h(lshell)* &
380 : orbtramat(lshell)%slm(is, ic)* &
381 2356816 : mpole_rho%Qlm_tot(iso)
382 :
383 : END DO
384 : END DO
385 : END DO ! lshell
386 : END DO ! iat
387 : END IF
388 : !MI Get rid of full gapw
389 :
390 36308 : CALL timestop(handle)
391 :
392 72616 : END SUBROUTINE calculate_rho0_atom
393 :
394 : ! **************************************************************************************************
395 : !> \brief ...
396 : !> \param local_rho_set ...
397 : !> \param qs_env ...
398 : !> \param gapw_control ...
399 : !> \param zcore ...
400 : ! **************************************************************************************************
401 6504 : SUBROUTINE init_rho0(local_rho_set, qs_env, gapw_control, zcore)
402 :
403 : TYPE(local_rho_type), POINTER :: local_rho_set
404 : TYPE(qs_environment_type), POINTER :: qs_env
405 : TYPE(gapw_control_type), POINTER :: gapw_control
406 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: zcore
407 :
408 : CHARACTER(len=*), PARAMETER :: routineN = 'init_rho0'
409 :
410 : CHARACTER(LEN=default_string_length) :: unit_str
411 : INTEGER :: handle, iat, iatom, ikind, l, l_rho1_max, l_rho1_max_nuc, laddg, lmaxg, maxl, &
412 : maxl_nuc, maxnset, maxso, maxso_nuc, nat, natom, nchan_c, nchan_s, nkind, nr, nset, &
413 : nset_nuc, nsotot, nsotot_nuc, output_unit
414 2168 : INTEGER, DIMENSION(:), POINTER :: atom_list
415 : LOGICAL :: cneo_potential_present, paw_atom
416 : REAL(KIND=dp) :: alpha_core, eps_Vrho0, max_rpgf0_s, radius, rc_min, rc_orb, &
417 : total_rho_core_rspace, total_rho_nuc_cneo_rspace, zeff
418 2168 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
419 : TYPE(cneo_potential_type), POINTER :: cneo_potential
420 : TYPE(cp_logger_type), POINTER :: logger
421 : TYPE(grid_atom_type), POINTER :: grid_atom
422 : TYPE(gto_basis_set_type), POINTER :: basis_1c, nuc_basis, nuc_soft_basis
423 : TYPE(harmonics_atom_type), POINTER :: harmonics
424 2168 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
425 2168 : TYPE(rho0_atom_type), DIMENSION(:), POINTER :: rho0_atom_set
426 : TYPE(rho0_mpole_type), POINTER :: rho0_mpole
427 2168 : TYPE(rhoz_cneo_type), DIMENSION(:), POINTER :: rhoz_cneo_set
428 2168 : TYPE(rhoz_type), DIMENSION(:), POINTER :: rhoz_set
429 : TYPE(section_vals_type), POINTER :: dft_section
430 :
431 2168 : CALL timeset(routineN, handle)
432 :
433 2168 : NULLIFY (logger)
434 2168 : logger => cp_get_default_logger()
435 :
436 2168 : NULLIFY (qs_kind_set)
437 2168 : NULLIFY (atomic_kind_set)
438 2168 : NULLIFY (harmonics)
439 2168 : NULLIFY (basis_1c)
440 2168 : NULLIFY (rho0_mpole)
441 2168 : NULLIFY (rho0_atom_set)
442 2168 : NULLIFY (rhoz_set)
443 2168 : NULLIFY (cneo_potential)
444 2168 : NULLIFY (nuc_basis)
445 2168 : NULLIFY (nuc_soft_basis)
446 2168 : NULLIFY (rhoz_cneo_set)
447 :
448 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
449 2168 : atomic_kind_set=atomic_kind_set)
450 :
451 2168 : nkind = SIZE(atomic_kind_set)
452 2168 : eps_Vrho0 = gapw_control%eps_Vrho0
453 :
454 : ! Initialize rhoz total to zero
455 : ! in gapw rhoz is calculated on local the lebedev grids
456 2168 : total_rho_core_rspace = 0.0_dp
457 2168 : total_rho_nuc_cneo_rspace = 0.0_dp
458 :
459 2168 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
460 :
461 : ! Initialize the multipole and the compensation charge type
462 2168 : CALL allocate_rho0_mpole(rho0_mpole)
463 2168 : CALL allocate_rho0_atom(rho0_atom_set, natom)
464 :
465 : ! Allocate the multipole set
466 2168 : CALL allocate_multipoles(rho0_mpole%mp_rho, natom, rho0_mpole%mp_gau, nkind)
467 :
468 : ! Allocate the core density on the radial grid for each kind: rhoz_set
469 2168 : CALL allocate_rhoz(rhoz_set, nkind)
470 :
471 : ! For each kind, determine the max l for the compensation charge density
472 2168 : lmaxg = gapw_control%lmax_rho0
473 2168 : laddg = gapw_control%ladd_rho0
474 :
475 2168 : CALL reallocate(rho0_mpole%lmax0_kind, 1, nkind)
476 :
477 2168 : rho0_mpole%lmax_0 = 0
478 2168 : rc_min = 100.0_dp
479 2168 : maxnset = 0
480 6448 : DO ikind = 1, nkind
481 4280 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
482 : CALL get_qs_kind(qs_kind_set(ikind), &
483 : ngrid_rad=nr, &
484 : grid_atom=grid_atom, &
485 : harmonics=harmonics, &
486 : paw_atom=paw_atom, &
487 : hard0_radius=rc_orb, &
488 : zeff=zeff, &
489 4280 : cneo_potential=cneo_potential)
490 : CALL get_qs_kind(qs_kind_set(ikind), &
491 4280 : basis_set=basis_1c, basis_type="GAPW_1C")
492 :
493 4280 : IF (ASSOCIATED(cneo_potential)) THEN
494 8 : IF (PRESENT(zcore)) THEN
495 0 : IF (zcore == 0.0_dp) THEN
496 0 : CPABORT("Electronic TDDFT with CNEO quantum nuclei is not implemented.")
497 : END IF
498 : END IF
499 8 : CPASSERT(paw_atom)
500 8 : NULLIFY (nuc_basis, nuc_soft_basis)
501 : CALL get_qs_kind(qs_kind_set(ikind), &
502 8 : basis_set=nuc_basis, basis_type="NUC")
503 : CALL get_qs_kind(qs_kind_set(ikind), &
504 8 : basis_set=nuc_soft_basis, basis_type="NUC_SOFT")
505 8 : alpha_core = 1.0_dp
506 : ELSE
507 4272 : CALL get_qs_kind(qs_kind_set(ikind), alpha_core_charge=alpha_core)
508 : END IF
509 :
510 : ! Set charge distribution of ionic cores to zero when computing the response-density
511 4280 : IF (PRESENT(zcore)) zeff = zcore
512 :
513 : CALL get_gto_basis_set(gto_basis_set=basis_1c, &
514 : maxl=maxl, &
515 4280 : maxso=maxso, nset=nset)
516 :
517 4280 : maxnset = MAX(maxnset, nset)
518 :
519 4280 : l_rho1_max = indso(1, harmonics%max_iso_not0)
520 :
521 4280 : maxl_nuc = -1
522 4280 : maxso_nuc = 0
523 4280 : nset_nuc = 0
524 4280 : l_rho1_max_nuc = -1
525 4280 : IF (ASSOCIATED(cneo_potential)) THEN
526 : CALL get_gto_basis_set(gto_basis_set=nuc_basis, &
527 : maxl=maxl_nuc, &
528 8 : maxso=maxso_nuc, nset=nset_nuc)
529 : ! Initialize CNEO potential internals
530 : CALL init_cneo_potential_internals(cneo_potential, nuc_basis, nuc_soft_basis, &
531 8 : gapw_control, grid_atom)
532 8 : l_rho1_max_nuc = indso(1, cneo_potential%harmonics%max_iso_not0)
533 : END IF
534 :
535 4280 : IF (paw_atom) THEN
536 : rho0_mpole%lmax0_kind(ikind) = MIN(2*MAX(maxl, maxl_nuc), &
537 : MAX(l_rho1_max, l_rho1_max_nuc), &
538 3926 : MAX(maxl, maxl_nuc) + laddg, lmaxg)
539 : ELSE
540 354 : rho0_mpole%lmax0_kind(ikind) = 0
541 : END IF
542 :
543 4280 : CALL set_qs_kind(qs_kind_set(ikind), lmax_rho0=rho0_mpole%lmax0_kind(ikind))
544 :
545 4280 : rho0_mpole%lmax_0 = MAX(rho0_mpole%lmax_0, rho0_mpole%lmax0_kind(ikind))
546 4280 : rc_min = MIN(rc_min, rc_orb)
547 :
548 4280 : nchan_s = nsoset(rho0_mpole%lmax0_kind(ikind))
549 4280 : nchan_c = ncoset(rho0_mpole%lmax0_kind(ikind))
550 4280 : nsotot = maxso*nset
551 4280 : nsotot_nuc = maxso_nuc*nset_nuc
552 :
553 11328 : DO iat = 1, nat
554 7048 : iatom = atom_list(iat)
555 : ! Allocate the multipole for rho1_h rho1_s and rho_z
556 7048 : CALL initialize_mpole_rho(rho0_mpole%mp_rho(iatom), nchan_s, nchan_c, zeff)
557 : ! Allocate the radial part of rho0_h and rho0_s
558 : ! This is calculated on the radial grid centered at the atomic position
559 11328 : CALL allocate_rho0_atom_rad(rho0_atom_set(iatom), nr, nchan_s)
560 : END DO
561 :
562 4280 : IF (paw_atom) THEN
563 : ! Calculate multipoles given by the product of 2 primitives Qlm_gg
564 : CALL calculate_mpole_gau(rho0_mpole%mp_gau(ikind)%Qlm_gg, &
565 3926 : basis_1c, harmonics, nchan_s, nsotot)
566 : END IF
567 :
568 15008 : IF (ASSOCIATED(cneo_potential)) THEN
569 8 : rho0_mpole%do_cneo = .TRUE.
570 : ! Calculate multipoles given by the product of two nuclear primitives Qlm_gg
571 : CALL calculate_mpole_gau(cneo_potential%Qlm_gg, nuc_basis, &
572 8 : cneo_potential%harmonics, nchan_s, nsotot_nuc)
573 : ! initial CNEO quantum nuclear charge density is a simple Zeff sum,
574 : ! but it will be calculated from numerical integration during SCF
575 8 : total_rho_nuc_cneo_rspace = total_rho_nuc_cneo_rspace - zeff*nat
576 : ELSE
577 : ! Calculate the core density rhoz
578 : ! exp(-alpha_c**2 r**2)Z(alpha_c**2/pi)**(3/2)
579 : ! on the logarithmic radial grid
580 : ! WARNING: alpha_core_charge = alpha_c**2
581 : CALL calculate_rhoz(rhoz_set(ikind), grid_atom, alpha_core, zeff, &
582 4272 : nat, total_rho_core_rspace, harmonics)
583 : END IF
584 : END DO ! ikind
585 2168 : total_rho_core_rspace = -total_rho_core_rspace
586 2168 : total_rho_nuc_cneo_rspace = -total_rho_nuc_cneo_rspace
587 :
588 : ! Allocate internals for quantum nuclear densities, if requested
589 2168 : CALL get_qs_kind_set(qs_kind_set, cneo_potential_present=cneo_potential_present)
590 2168 : IF (cneo_potential_present) THEN
591 : CALL allocate_rhoz_cneo_internals(rhoz_cneo_set, atomic_kind_set, &
592 8 : qs_kind_set, qs_env)
593 : END IF
594 :
595 2168 : IF (gapw_control%alpha0_hard_from_input) THEN
596 : ! The exponent for the compensation charge rho0_hard is read from input
597 110 : rho0_mpole%zet0_h = gapw_control%alpha0_hard
598 : ELSE
599 : ! Calculate the exponent for the compensation charge rho0_hard
600 2058 : rho0_mpole%zet0_h = 0.1_dp
601 184222 : DO
602 186280 : radius = exp_radius(rho0_mpole%lmax_0, rho0_mpole%zet0_h, eps_Vrho0, 1.0_dp)
603 186280 : IF (radius <= rc_min) EXIT
604 184222 : rho0_mpole%zet0_h = rho0_mpole%zet0_h + 0.1_dp
605 : END DO
606 :
607 : END IF
608 :
609 : ! Allocate and calculate the normalization factors for g0_lm_h and g0_lm_s
610 2168 : CALL reallocate(rho0_mpole%norm_g0l_h, 0, rho0_mpole%lmax_0)
611 8664 : DO l = 0, rho0_mpole%lmax_0
612 : rho0_mpole%norm_g0l_h(l) = (2._dp*l + 1._dp)/ &
613 8664 : (fourpi*gaussint_sph(rho0_mpole%zet0_h, 2*l))
614 : END DO
615 :
616 : ! Allocate and Initialize the g0 gaussians used to build the compensation density
617 : ! and calculate the interaction radii
618 2168 : max_rpgf0_s = 0.0_dp
619 6448 : DO ikind = 1, nkind
620 4280 : CALL get_qs_kind(qs_kind_set(ikind), grid_atom=grid_atom)
621 4280 : CALL calculate_g0(rho0_mpole, grid_atom, ikind)
622 6448 : CALL interaction_radii_g0(rho0_mpole, ikind, eps_Vrho0, max_rpgf0_s)
623 : END DO
624 2168 : rho0_mpole%max_rpgf0_s = max_rpgf0_s
625 :
626 : CALL set_local_rho(local_rho_set, rho0_atom_set=rho0_atom_set, rho0_mpole=rho0_mpole, &
627 2168 : rhoz_set=rhoz_set, rhoz_cneo_set=rhoz_cneo_set)
628 2168 : local_rho_set%rhoz_tot = total_rho_core_rspace
629 2168 : local_rho_set%rhoz_cneo_tot = total_rho_nuc_cneo_rspace
630 :
631 2168 : dft_section => section_vals_get_subs_vals(qs_env%input, "DFT")
632 : output_unit = cp_print_key_unit_nr(logger, dft_section, "PRINT%GAPW%RHO0_INFORMATION", &
633 2168 : extension=".Log")
634 2168 : CALL section_vals_val_get(dft_section, "PRINT%GAPW%RHO0_INFORMATION%UNIT", c_val=unit_str)
635 2168 : IF (output_unit > 0) THEN
636 1 : CALL write_rho0_info(rho0_mpole, unit_str, output_unit)
637 : END IF
638 : CALL cp_print_key_finished_output(output_unit, logger, dft_section, &
639 2168 : "PRINT%GAPW%RHO0_INFORMATION")
640 :
641 2168 : CALL timestop(handle)
642 :
643 2168 : END SUBROUTINE init_rho0
644 :
645 : ! **************************************************************************************************
646 : !> \brief ...
647 : !> \param rho0_mpole ...
648 : !> \param ik ...
649 : !> \param eps_Vrho0 ...
650 : !> \param max_rpgf0_s ...
651 : ! **************************************************************************************************
652 4280 : SUBROUTINE interaction_radii_g0(rho0_mpole, ik, eps_Vrho0, max_rpgf0_s)
653 :
654 : TYPE(rho0_mpole_type), POINTER :: rho0_mpole
655 : INTEGER, INTENT(IN) :: ik
656 : REAL(KIND=dp), INTENT(IN) :: eps_Vrho0
657 : REAL(KIND=dp), INTENT(INOUT) :: max_rpgf0_s
658 :
659 : INTEGER :: l, lmax
660 : REAL(KIND=dp) :: r_h, z0_h
661 4280 : REAL(KIND=dp), DIMENSION(:), POINTER :: ng0_h
662 :
663 : CALL get_rho0_mpole(rho0_mpole, ikind=ik, l0_ikind=lmax, &
664 4280 : zet0_h=z0_h, norm_g0l_h=ng0_h)
665 4280 : r_h = 0.0_dp
666 16004 : DO l = 0, lmax
667 16004 : r_h = MAX(r_h, exp_radius(l, z0_h, eps_Vrho0, ng0_h(l), rlow=r_h))
668 : END DO
669 :
670 4280 : rho0_mpole%mp_gau(ik)%rpgf0_h = r_h
671 4280 : rho0_mpole%mp_gau(ik)%rpgf0_s = r_h
672 4280 : max_rpgf0_s = MAX(max_rpgf0_s, r_h)
673 :
674 4280 : END SUBROUTINE interaction_radii_g0
675 :
676 : END MODULE qs_rho0_methods
|
