Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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 6714 : 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 6714 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf
103 : REAL(KIND=dp) :: zet1, zet2
104 6714 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zet
105 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: my_CG
106 :
107 6714 : CALL timeset(routineN, handle)
108 :
109 6714 : NULLIFY (lmax, lmin, npgf, my_CG, zet)
110 :
111 : CALL reallocate(Qlm_gg, 1, nsotot, 1, nsotot, 1, &
112 6714 : 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 6714 : npgf=npgf, nset=nset, zet=zet, maxl=maxl)
117 :
118 6714 : max_s_harm = harmonics%max_s_harm
119 6714 : llmax = harmonics%llmax
120 :
121 40284 : ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm))
122 :
123 6714 : my_CG => harmonics%my_CG
124 :
125 6714 : m1 = 0
126 24872 : DO iset1 = 1, nset
127 : m2 = 0
128 82144 : DO iset2 = 1, nset
129 :
130 : CALL get_none0_cg_list(my_CG, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
131 63986 : max_s_harm, llmax, cg_list, cg_n_list, max_iso_not0_local)
132 :
133 63986 : n1 = nsoset(lmax(iset1))
134 219238 : DO ipgf1 = 1, npgf(iset1)
135 155252 : zet1 = zet(ipgf1, iset1)
136 :
137 155252 : n2 = nsoset(lmax(iset2))
138 635028 : DO ipgf2 = 1, npgf(iset2)
139 415790 : zet2 = zet(ipgf2, iset2)
140 :
141 3113504 : DO iso = 1, MIN(nchannels, max_iso_not0_local)
142 2542462 : l = indso(1, iso)
143 7922482 : DO icg = 1, cg_n_list(iso)
144 4964230 : iso1 = cg_list(1, icg, iso)
145 4964230 : iso2 = cg_list(2, icg, iso)
146 :
147 4964230 : l1 = indso(1, iso1)
148 4964230 : l2 = indso(1, iso2)
149 4964230 : ig1 = iso1 + n1*(ipgf1 - 1) + m1
150 4964230 : ig2 = iso2 + n2*(ipgf2 - 1) + m2
151 :
152 : Qlm_gg(ig1, ig2, iso) = fourpi/(2._dp*l + 1._dp)* &
153 7506692 : 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 146130 : m2 = m2 + maxso
160 : END DO ! iset2
161 24872 : m1 = m1 + maxso
162 : END DO ! iset1
163 :
164 6714 : DEALLOCATE (cg_list, cg_n_list)
165 :
166 6714 : CALL timestop(handle)
167 13428 : 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 53718 : SUBROUTINE calculate_rho0_atom(gapw_control, rho_atom_set, rhoz_cneo_set, rho0_atom_set, &
183 53718 : 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 53718 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cpc_h_nuc, cpc_s_nuc
203 53718 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: cpc_ah, cpc_as
204 53718 : REAL(KIND=dp), DIMENSION(:), POINTER :: norm_g0l_h
205 53718 : 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 53718 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: cpc_h, cpc_s
212 : TYPE(rho_atom_type), POINTER :: rho_atom
213 :
214 53718 : CALL timeset(routineN, handle)
215 :
216 53718 : NULLIFY (mpole_gau)
217 53718 : NULLIFY (mpole_rho)
218 53718 : NULLIFY (g0_h, vg0_h, g_atom)
219 53718 : NULLIFY (norm_g0l_h, harmonics)
220 53718 : 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 53718 : 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 53718 : cneo_potential=cneo_potential)
230 :
231 53718 : nr = g_atom%nr
232 :
233 : ! Set density coefficient to zero before the calculation
234 140256 : DO iat = 1, natom
235 86538 : iatom = a_list(iat)
236 34907152 : rho0_atom_set(iatom)%rho0_rad_h%r_coef = 0.0_dp
237 737572 : rho0_mp%mp_rho(iatom)%Qlm_tot = 0.0_dp
238 : ! When no nuclear density matrix is available, use spherical zeff
239 86538 : IF (.NOT. ASSOCIATED(cneo_potential)) THEN
240 86446 : 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 259614 : rho0_mp%mp_rho(iatom)%Q0 = 0.0_dp
247 871022 : rho0_mp%mp_rho(iatom)%Qlm_car = 0.0_dp
248 : END DO
249 :
250 53718 : IF (.NOT. (.NOT. paw_atom .AND. gapw_control%nopaw_as_gpw)) THEN
251 133152 : DO iat = 1, natom
252 81376 : iatom = a_list(iat)
253 81376 : mpole_rho => rho0_mp%mp_rho(iatom)
254 81376 : rho_atom => rho_atom_set(iatom)
255 :
256 81376 : IF (paw_atom) THEN
257 81376 : NULLIFY (cpc_h, cpc_s)
258 81376 : CALL get_rho_atom(rho_atom=rho_atom, cpc_h=cpc_h, cpc_s=cpc_s)
259 81376 : nspins = SIZE(cpc_h)
260 81376 : nsotot = SIZE(mpole_gau%Qlm_gg, 1)
261 81376 : IF (ALLOCATED(cpc_ah)) THEN
262 29600 : IF (SIZE(cpc_ah, 1) /= nsotot .OR. SIZE(cpc_ah, 2) /= nsotot .OR. &
263 0 : SIZE(cpc_ah, 3) /= nspins) DEALLOCATE (cpc_ah)
264 : END IF
265 288480 : IF (.NOT. ALLOCATED(cpc_ah)) ALLOCATE (cpc_ah(nsotot, nsotot, nspins))
266 303316208 : cpc_ah = 0._dp
267 81376 : IF (ALLOCATED(cpc_as)) THEN
268 29600 : IF (SIZE(cpc_as, 1) /= nsotot .OR. SIZE(cpc_as, 2) /= nsotot .OR. &
269 0 : SIZE(cpc_as, 3) /= nspins) DEALLOCATE (cpc_as)
270 : END IF
271 288480 : IF (.NOT. ALLOCATED(cpc_as)) ALLOCATE (cpc_as(nsotot, nsotot, nspins))
272 303316208 : cpc_as = 0._dp
273 172908 : DO ispin = 1, nspins
274 91532 : CALL prj_scatter(cpc_h(ispin)%r_coef, cpc_ah(:, :, ispin), qs_kind)
275 172908 : CALL prj_scatter(cpc_s(ispin)%r_coef, cpc_as(:, :, ispin), qs_kind)
276 : END DO
277 81376 : nsotot_nuc = 0
278 81376 : IF (ASSOCIATED(cneo_potential)) THEN
279 92 : IF (rhoz_cneo_set(iatom)%ready) THEN
280 78 : nsotot_nuc = SIZE(cneo_potential%Qlm_gg, 1)
281 78 : IF (ALLOCATED(cpc_h_nuc)) THEN
282 38 : IF (SIZE(cpc_h_nuc, 1) /= nsotot_nuc .OR. &
283 0 : SIZE(cpc_h_nuc, 2) /= nsotot_nuc) DEALLOCATE (cpc_h_nuc)
284 : END IF
285 198 : IF (.NOT. ALLOCATED(cpc_h_nuc)) ALLOCATE (cpc_h_nuc(nsotot_nuc, nsotot_nuc))
286 78 : IF (ALLOCATED(cpc_s_nuc)) THEN
287 38 : IF (SIZE(cpc_s_nuc, 1) /= nsotot_nuc .OR. &
288 0 : SIZE(cpc_s_nuc, 2) /= nsotot_nuc) DEALLOCATE (cpc_s_nuc)
289 : END IF
290 198 : IF (.NOT. ALLOCATED(cpc_s_nuc)) ALLOCATE (cpc_s_nuc(nsotot_nuc, nsotot_nuc))
291 518154 : cpc_h_nuc = 0._dp
292 518154 : cpc_s_nuc = 0._dp
293 : CALL cneo_scatter(rhoz_cneo_set(iatom)%cpc_h, cpc_h_nuc, cneo_potential%npsgf, &
294 78 : cneo_potential%n2oindex)
295 : CALL cneo_scatter(rhoz_cneo_set(iatom)%cpc_s, cpc_s_nuc, cneo_potential%npsgf, &
296 78 : cneo_potential%n2oindex)
297 : END IF
298 : END IF
299 :
300 : ! Total charge (hard-soft) at atom
301 : IF (paw_atom) THEN
302 172908 : DO ispin = 1, nspins
303 : mpole_rho%Q0(ispin) = (trace_r_AxB(mpole_gau%Qlm_gg(:, :, 1), nsotot, &
304 : cpc_ah(:, :, ispin), nsotot, nsotot, nsotot) &
305 : - trace_r_AxB(mpole_gau%Qlm_gg(:, :, 1), nsotot, &
306 172908 : cpc_as(:, :, ispin), nsotot, nsotot, nsotot))/SQRT(fourpi)
307 : END DO
308 : END IF
309 : ! Multipoles of local charge distribution
310 727248 : DO iso = 1, MIN(nsoset(lmax0), harmonics%max_iso_not0)
311 81376 : IF (paw_atom) THEN
312 645872 : mpole_rho%Qlm_h(iso) = 0.0_dp
313 645872 : mpole_rho%Qlm_s(iso) = 0.0_dp
314 :
315 1363500 : DO ispin = 1, nspins
316 : mpole_rho%Qlm_h(iso) = mpole_rho%Qlm_h(iso) + &
317 : trace_r_AxB(mpole_gau%Qlm_gg(:, :, iso), nsotot, &
318 717628 : cpc_ah(:, :, ispin), nsotot, nsotot, nsotot)
319 : mpole_rho%Qlm_s(iso) = mpole_rho%Qlm_s(iso) + &
320 : trace_r_AxB(mpole_gau%Qlm_gg(:, :, iso), nsotot, &
321 1363500 : cpc_as(:, :, ispin), nsotot, nsotot, nsotot)
322 : END DO ! ispin
323 :
324 : mpole_rho%Qlm_tot(iso) = mpole_rho%Qlm_tot(iso) + &
325 645872 : mpole_rho%Qlm_h(iso) - mpole_rho%Qlm_s(iso)
326 : END IF
327 : END DO ! iso
328 :
329 : ! Multipoles of CNEO quantum nuclear charge distribuition
330 81376 : IF (ASSOCIATED(cneo_potential)) THEN
331 92 : IF (rhoz_cneo_set(iatom)%ready) THEN
332 780 : DO iso = 1, MIN(nsoset(lmax0), cneo_potential%harmonics%max_iso_not0)
333 : mpole_rho%Qlm_tot(iso) = mpole_rho%Qlm_tot(iso) - cneo_potential%zeff* &
334 : trace_r_AxB(cneo_potential%Qlm_gg(:, :, iso), nsotot_nuc, &
335 : cpc_h_nuc - cpc_s_nuc, nsotot_nuc, &
336 4663464 : nsotot_nuc, nsotot_nuc)
337 : END DO ! iso
338 : END IF
339 : END IF
340 : END IF
341 :
342 780966 : DO iso = 1, nsoset(lmax0)
343 645872 : l = indso(1, iso)
344 : rho0_atom_set(iatom)%rho0_rad_h%r_coef(1:nr, iso) = &
345 68469312 : g0_h(1:nr, l)*mpole_rho%Qlm_tot(iso)
346 : rho0_atom_set(iatom)%vrho0_rad_h%r_coef(1:nr, iso) = &
347 68469312 : vg0_h(1:nr, l)*mpole_rho%Qlm_tot(iso)
348 :
349 : ! When CNEO is enabled, it is possible for rho0 to have a higher angualr momentum
350 : ! than that of the electronic density. In that case, the nuclear density must have
351 : ! a higher angular momentum, but cneo_potential%harmonics%slm_int is not initialized.
352 : ! For higher angular momenta, simply use the fact that slm_int(iso>1)=0
353 727248 : IF (iso <= harmonics%max_iso_not0) THEN
354 : sum1 = 0.0_dp
355 34557592 : DO ir = 1, nr
356 : sum1 = sum1 + g_atom%wr(ir)* &
357 34557592 : rho0_atom_set(iatom)%rho0_rad_h%r_coef(ir, iso)
358 : END DO
359 645872 : rho0_h_tot = rho0_h_tot + sum1*harmonics%slm_int(iso)
360 : END IF
361 : END DO ! iso
362 : END DO ! iat
363 : END IF
364 :
365 53718 : IF (ALLOCATED(cpc_ah)) DEALLOCATE (cpc_ah)
366 53718 : IF (ALLOCATED(cpc_as)) DEALLOCATE (cpc_as)
367 53718 : IF (ALLOCATED(cpc_h_nuc)) DEALLOCATE (cpc_h_nuc)
368 53718 : IF (ALLOCATED(cpc_s_nuc)) DEALLOCATE (cpc_s_nuc)
369 :
370 : ! Transform the coefficinets from spherical to Cartesian
371 53718 : IF (.NOT. paw_atom .AND. gapw_control%nopaw_as_gpw) THEN
372 7104 : DO iat = 1, natom
373 5162 : iatom = a_list(iat)
374 5162 : mpole_rho => rho0_mp%mp_rho(iatom)
375 :
376 12266 : DO lshell = 0, lmax0
377 15486 : DO ic = 1, nco(lshell)
378 5162 : ico = ic + ncoset(lshell - 1)
379 10324 : mpole_rho%Qlm_car(ico) = 0.0_dp
380 : END DO
381 : END DO
382 : END DO
383 : ELSE
384 133152 : DO iat = 1, natom
385 81376 : iatom = a_list(iat)
386 81376 : mpole_rho => rho0_mp%mp_rho(iatom)
387 353032 : DO lshell = 0, lmax0
388 1026860 : DO ic = 1, nco(lshell)
389 725604 : ico = ic + ncoset(lshell - 1)
390 725604 : mpole_rho%Qlm_car(ico) = 0.0_dp
391 725604 : lx = indco(1, ico)
392 725604 : ly = indco(2, ico)
393 725604 : lz = indco(3, ico)
394 3945148 : DO is = 1, nso(lshell)
395 2999664 : iso = is + nsoset(lshell - 1)
396 : mpole_rho%Qlm_car(ico) = mpole_rho%Qlm_car(ico) + &
397 : norm_g0l_h(lshell)* &
398 : orbtramat(lshell)%slm(is, ic)* &
399 3725268 : mpole_rho%Qlm_tot(iso)
400 :
401 : END DO
402 : END DO
403 : END DO ! lshell
404 : END DO ! iat
405 : END IF
406 : !MI Get rid of full gapw
407 :
408 53718 : CALL timestop(handle)
409 :
410 107436 : END SUBROUTINE calculate_rho0_atom
411 :
412 : ! **************************************************************************************************
413 : !> \brief ...
414 : !> \param local_rho_set ...
415 : !> \param qs_env ...
416 : !> \param gapw_control ...
417 : !> \param zcore ...
418 : ! **************************************************************************************************
419 10614 : SUBROUTINE init_rho0(local_rho_set, qs_env, gapw_control, zcore)
420 :
421 : TYPE(local_rho_type), POINTER :: local_rho_set
422 : TYPE(qs_environment_type), POINTER :: qs_env
423 : TYPE(gapw_control_type), POINTER :: gapw_control
424 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: zcore
425 :
426 : CHARACTER(len=*), PARAMETER :: routineN = 'init_rho0'
427 :
428 : CHARACTER(LEN=default_string_length) :: unit_str
429 : INTEGER :: handle, iat, iatom, ikind, l, l_rho1_max, l_rho1_max_nuc, laddg, lmaxg, maxl, &
430 : maxl_nuc, maxnset, maxso, maxso_nuc, nat, natom, nchan_c, nchan_s, nkind, nr, nset, &
431 : nset_nuc, nsotot, nsotot_nuc, output_unit
432 3538 : INTEGER, DIMENSION(:), POINTER :: atom_list
433 : LOGICAL :: cneo_potential_present, paw_atom
434 : REAL(KIND=dp) :: alpha_core, eps_Vrho0, max_rpgf0_s, radius, rc_min, rc_orb, &
435 : total_rho_core_rspace, total_rho_nuc_cneo_rspace, zeff
436 3538 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
437 : TYPE(cneo_potential_type), POINTER :: cneo_potential
438 : TYPE(cp_logger_type), POINTER :: logger
439 : TYPE(grid_atom_type), POINTER :: grid_atom
440 : TYPE(gto_basis_set_type), POINTER :: basis_1c, nuc_basis, nuc_soft_basis
441 : TYPE(harmonics_atom_type), POINTER :: harmonics
442 3538 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
443 3538 : TYPE(rho0_atom_type), DIMENSION(:), POINTER :: rho0_atom_set
444 : TYPE(rho0_mpole_type), POINTER :: rho0_mpole
445 3538 : TYPE(rhoz_cneo_type), DIMENSION(:), POINTER :: rhoz_cneo_set
446 3538 : TYPE(rhoz_type), DIMENSION(:), POINTER :: rhoz_set
447 : TYPE(section_vals_type), POINTER :: dft_section
448 :
449 3538 : CALL timeset(routineN, handle)
450 :
451 3538 : NULLIFY (logger)
452 3538 : logger => cp_get_default_logger()
453 :
454 3538 : NULLIFY (qs_kind_set)
455 3538 : NULLIFY (atomic_kind_set)
456 3538 : NULLIFY (harmonics)
457 3538 : NULLIFY (basis_1c)
458 3538 : NULLIFY (rho0_mpole)
459 3538 : NULLIFY (rho0_atom_set)
460 3538 : NULLIFY (rhoz_set)
461 3538 : NULLIFY (cneo_potential)
462 3538 : NULLIFY (nuc_basis)
463 3538 : NULLIFY (nuc_soft_basis)
464 3538 : NULLIFY (rhoz_cneo_set)
465 :
466 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
467 3538 : atomic_kind_set=atomic_kind_set)
468 :
469 3538 : nkind = SIZE(atomic_kind_set)
470 3538 : eps_Vrho0 = gapw_control%eps_Vrho0
471 :
472 : ! Initialize rhoz total to zero
473 : ! in gapw rhoz is calculated on local the lebedev grids
474 3538 : total_rho_core_rspace = 0.0_dp
475 3538 : total_rho_nuc_cneo_rspace = 0.0_dp
476 :
477 3538 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
478 :
479 : ! Initialize the multipole and the compensation charge type
480 3538 : CALL allocate_rho0_mpole(rho0_mpole)
481 3538 : CALL allocate_rho0_atom(rho0_atom_set, natom)
482 :
483 : ! Allocate the multipole set
484 3538 : CALL allocate_multipoles(rho0_mpole%mp_rho, natom, rho0_mpole%mp_gau, nkind)
485 :
486 : ! Allocate the core density on the radial grid for each kind: rhoz_set
487 3538 : CALL allocate_rhoz(rhoz_set, nkind)
488 :
489 : ! For each kind, determine the max l for the compensation charge density
490 3538 : lmaxg = gapw_control%lmax_rho0
491 3538 : laddg = gapw_control%ladd_rho0
492 :
493 3538 : CALL reallocate(rho0_mpole%lmax0_kind, 1, nkind)
494 :
495 3538 : rho0_mpole%lmax_0 = 0
496 3538 : rc_min = 100.0_dp
497 3538 : maxnset = 0
498 10662 : DO ikind = 1, nkind
499 7124 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
500 : CALL get_qs_kind(qs_kind_set(ikind), &
501 : ngrid_rad=nr, &
502 : grid_atom=grid_atom, &
503 : harmonics=harmonics, &
504 : paw_atom=paw_atom, &
505 : hard0_radius=rc_orb, &
506 : zeff=zeff, &
507 7124 : cneo_potential=cneo_potential)
508 : CALL get_qs_kind(qs_kind_set(ikind), &
509 7124 : basis_set=basis_1c, basis_type="GAPW_1C")
510 :
511 7124 : IF (ASSOCIATED(cneo_potential)) THEN
512 8 : IF (PRESENT(zcore)) THEN
513 0 : IF (zcore == 0.0_dp) THEN
514 0 : CPABORT("Electronic TDDFT with CNEO quantum nuclei is not implemented.")
515 : END IF
516 : END IF
517 8 : CPASSERT(paw_atom)
518 8 : NULLIFY (nuc_basis, nuc_soft_basis)
519 : CALL get_qs_kind(qs_kind_set(ikind), &
520 8 : basis_set=nuc_basis, basis_type="NUC")
521 : CALL get_qs_kind(qs_kind_set(ikind), &
522 8 : basis_set=nuc_soft_basis, basis_type="NUC_SOFT")
523 8 : alpha_core = 1.0_dp
524 : ELSE
525 7116 : CALL get_qs_kind(qs_kind_set(ikind), alpha_core_charge=alpha_core)
526 : END IF
527 :
528 : ! Set charge distribution of ionic cores to zero when computing the response-density
529 7124 : IF (PRESENT(zcore)) zeff = zcore
530 :
531 : CALL get_gto_basis_set(gto_basis_set=basis_1c, &
532 : maxl=maxl, &
533 7124 : maxso=maxso, nset=nset)
534 :
535 7124 : maxnset = MAX(maxnset, nset)
536 :
537 7124 : l_rho1_max = indso(1, harmonics%max_iso_not0)
538 :
539 7124 : maxl_nuc = -1
540 7124 : maxso_nuc = 0
541 7124 : nset_nuc = 0
542 7124 : l_rho1_max_nuc = -1
543 7124 : IF (ASSOCIATED(cneo_potential)) THEN
544 : CALL get_gto_basis_set(gto_basis_set=nuc_basis, &
545 : maxl=maxl_nuc, &
546 8 : maxso=maxso_nuc, nset=nset_nuc)
547 : ! Initialize CNEO potential internals
548 : CALL init_cneo_potential_internals(cneo_potential, nuc_basis, nuc_soft_basis, &
549 8 : gapw_control, grid_atom)
550 8 : l_rho1_max_nuc = indso(1, cneo_potential%harmonics%max_iso_not0)
551 : END IF
552 :
553 7124 : IF (paw_atom) THEN
554 : rho0_mpole%lmax0_kind(ikind) = MIN(2*MAX(maxl, maxl_nuc), &
555 : MAX(l_rho1_max, l_rho1_max_nuc), &
556 6706 : MAX(maxl, maxl_nuc) + laddg, lmaxg)
557 : ELSE
558 418 : rho0_mpole%lmax0_kind(ikind) = 0
559 : END IF
560 :
561 7124 : CALL set_qs_kind(qs_kind_set(ikind), lmax_rho0=rho0_mpole%lmax0_kind(ikind))
562 :
563 7124 : rho0_mpole%lmax_0 = MAX(rho0_mpole%lmax_0, rho0_mpole%lmax0_kind(ikind))
564 7124 : rc_min = MIN(rc_min, rc_orb)
565 :
566 7124 : nchan_s = nsoset(rho0_mpole%lmax0_kind(ikind))
567 7124 : nchan_c = ncoset(rho0_mpole%lmax0_kind(ikind))
568 7124 : nsotot = maxso*nset
569 7124 : nsotot_nuc = maxso_nuc*nset_nuc
570 :
571 18332 : DO iat = 1, nat
572 11208 : iatom = atom_list(iat)
573 : ! Allocate the multipole for rho1_h rho1_s and rho_z
574 11208 : CALL initialize_mpole_rho(rho0_mpole%mp_rho(iatom), nchan_s, nchan_c, zeff)
575 : ! Allocate the radial part of rho0_h and rho0_s
576 : ! This is calculated on the radial grid centered at the atomic position
577 18332 : CALL allocate_rho0_atom_rad(rho0_atom_set(iatom), nr, nchan_s)
578 : END DO
579 :
580 7124 : IF (paw_atom) THEN
581 : ! Calculate multipoles given by the product of 2 primitives Qlm_gg
582 : CALL calculate_mpole_gau(rho0_mpole%mp_gau(ikind)%Qlm_gg, &
583 6706 : basis_1c, harmonics, nchan_s, nsotot)
584 : END IF
585 :
586 24910 : IF (ASSOCIATED(cneo_potential)) THEN
587 8 : rho0_mpole%do_cneo = .TRUE.
588 : ! Calculate multipoles given by the product of two nuclear primitives Qlm_gg
589 : CALL calculate_mpole_gau(cneo_potential%Qlm_gg, nuc_basis, &
590 8 : cneo_potential%harmonics, nchan_s, nsotot_nuc)
591 : ! initial CNEO quantum nuclear charge density is a simple Zeff sum,
592 : ! but it will be calculated from numerical integration during SCF
593 8 : total_rho_nuc_cneo_rspace = total_rho_nuc_cneo_rspace - zeff*nat
594 : ELSE
595 : ! Calculate the core density rhoz
596 : ! exp(-alpha_c**2 r**2)Z(alpha_c**2/pi)**(3/2)
597 : ! on the logarithmic radial grid
598 : ! WARNING: alpha_core_charge = alpha_c**2
599 : CALL calculate_rhoz(rhoz_set(ikind), grid_atom, alpha_core, zeff, &
600 7116 : nat, total_rho_core_rspace, harmonics)
601 : END IF
602 : END DO ! ikind
603 3538 : total_rho_core_rspace = -total_rho_core_rspace
604 3538 : total_rho_nuc_cneo_rspace = -total_rho_nuc_cneo_rspace
605 :
606 : ! Allocate internals for quantum nuclear densities, if requested
607 3538 : CALL get_qs_kind_set(qs_kind_set, cneo_potential_present=cneo_potential_present)
608 3538 : IF (cneo_potential_present) THEN
609 : CALL allocate_rhoz_cneo_internals(rhoz_cneo_set, atomic_kind_set, &
610 8 : qs_kind_set, qs_env)
611 : END IF
612 :
613 3538 : IF (gapw_control%alpha0_hard_from_input) THEN
614 : ! The exponent for the compensation charge rho0_hard is read from input
615 116 : rho0_mpole%zet0_h = gapw_control%alpha0_hard
616 : ELSE
617 : ! Calculate the exponent for the compensation charge rho0_hard
618 3422 : rho0_mpole%zet0_h = 0.1_dp
619 313392 : DO
620 316814 : radius = exp_radius(rho0_mpole%lmax_0, rho0_mpole%zet0_h, eps_Vrho0, 1.0_dp)
621 316814 : IF (radius <= rc_min) EXIT
622 313392 : rho0_mpole%zet0_h = rho0_mpole%zet0_h + 0.1_dp
623 : END DO
624 :
625 : END IF
626 :
627 : ! Allocate and calculate the normalization factors for g0_lm_h and g0_lm_s
628 3538 : CALL reallocate(rho0_mpole%norm_g0l_h, 0, rho0_mpole%lmax_0)
629 14320 : DO l = 0, rho0_mpole%lmax_0
630 : rho0_mpole%norm_g0l_h(l) = (2._dp*l + 1._dp)/ &
631 14320 : (fourpi*gaussint_sph(rho0_mpole%zet0_h, 2*l))
632 : END DO
633 :
634 : ! Allocate and Initialize the g0 gaussians used to build the compensation density
635 : ! and calculate the interaction radii
636 3538 : max_rpgf0_s = 0.0_dp
637 10662 : DO ikind = 1, nkind
638 7124 : CALL get_qs_kind(qs_kind_set(ikind), grid_atom=grid_atom)
639 7124 : CALL calculate_g0(rho0_mpole, grid_atom, ikind)
640 10662 : CALL interaction_radii_g0(rho0_mpole, ikind, eps_Vrho0, max_rpgf0_s)
641 : END DO
642 3538 : rho0_mpole%max_rpgf0_s = max_rpgf0_s
643 :
644 : CALL set_local_rho(local_rho_set, rho0_atom_set=rho0_atom_set, rho0_mpole=rho0_mpole, &
645 3538 : rhoz_set=rhoz_set, rhoz_cneo_set=rhoz_cneo_set)
646 3538 : local_rho_set%rhoz_tot = total_rho_core_rspace
647 3538 : local_rho_set%rhoz_cneo_tot = total_rho_nuc_cneo_rspace
648 :
649 3538 : dft_section => section_vals_get_subs_vals(qs_env%input, "DFT")
650 : output_unit = cp_print_key_unit_nr(logger, dft_section, "PRINT%GAPW%RHO0_INFORMATION", &
651 3538 : extension=".Log")
652 3538 : CALL section_vals_val_get(dft_section, "PRINT%GAPW%RHO0_INFORMATION%UNIT", c_val=unit_str)
653 3538 : IF (output_unit > 0) THEN
654 1 : CALL write_rho0_info(rho0_mpole, unit_str, output_unit)
655 : END IF
656 : CALL cp_print_key_finished_output(output_unit, logger, dft_section, &
657 3538 : "PRINT%GAPW%RHO0_INFORMATION")
658 :
659 3538 : CALL timestop(handle)
660 :
661 3538 : END SUBROUTINE init_rho0
662 :
663 : ! **************************************************************************************************
664 : !> \brief ...
665 : !> \param rho0_mpole ...
666 : !> \param ik ...
667 : !> \param eps_Vrho0 ...
668 : !> \param max_rpgf0_s ...
669 : ! **************************************************************************************************
670 7124 : SUBROUTINE interaction_radii_g0(rho0_mpole, ik, eps_Vrho0, max_rpgf0_s)
671 :
672 : TYPE(rho0_mpole_type), POINTER :: rho0_mpole
673 : INTEGER, INTENT(IN) :: ik
674 : REAL(KIND=dp), INTENT(IN) :: eps_Vrho0
675 : REAL(KIND=dp), INTENT(INOUT) :: max_rpgf0_s
676 :
677 : INTEGER :: l, lmax
678 : REAL(KIND=dp) :: r_h, z0_h
679 7124 : REAL(KIND=dp), DIMENSION(:), POINTER :: ng0_h
680 :
681 : CALL get_rho0_mpole(rho0_mpole, ikind=ik, l0_ikind=lmax, &
682 7124 : zet0_h=z0_h, norm_g0l_h=ng0_h)
683 7124 : r_h = 0.0_dp
684 27248 : DO l = 0, lmax
685 27248 : r_h = MAX(r_h, exp_radius(l, z0_h, eps_Vrho0, ng0_h(l), rlow=r_h))
686 : END DO
687 :
688 7124 : rho0_mpole%mp_gau(ik)%rpgf0_h = r_h
689 7124 : rho0_mpole%mp_gau(ik)%rpgf0_s = r_h
690 7124 : max_rpgf0_s = MAX(max_rpgf0_s, r_h)
691 :
692 7124 : END SUBROUTINE interaction_radii_g0
693 :
694 : END MODULE qs_rho0_methods
|
