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 : ! **************************************************************************************************
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_environment_types, ONLY: get_qs_env,&
39 : qs_environment_type
40 : USE qs_grid_atom, ONLY: grid_atom_type
41 : USE qs_harmonics_atom, ONLY: get_none0_cg_list,&
42 : harmonics_atom_type
43 : USE qs_kind_types, ONLY: get_qs_kind,&
44 : qs_kind_type,&
45 : set_qs_kind
46 : USE qs_local_rho_types, ONLY: allocate_rhoz,&
47 : calculate_rhoz,&
48 : local_rho_type,&
49 : rhoz_type,&
50 : set_local_rho
51 : USE qs_oce_methods, ONLY: prj_scatter
52 : USE qs_rho0_types, ONLY: &
53 : allocate_multipoles, allocate_rho0_atom, allocate_rho0_atom_rad, allocate_rho0_mpole, &
54 : calculate_g0, get_rho0_mpole, initialize_mpole_rho, mpole_gau_overlap, mpole_rho_atom, &
55 : rho0_atom_type, rho0_mpole_type, write_rho0_info
56 : USE qs_rho_atom_types, ONLY: get_rho_atom,&
57 : rho_atom_coeff,&
58 : rho_atom_type
59 : #include "./base/base_uses.f90"
60 :
61 : IMPLICIT NONE
62 :
63 : PRIVATE
64 :
65 : ! Global parameters (only in this module)
66 :
67 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_rho0_methods'
68 :
69 : ! Public subroutines
70 :
71 : PUBLIC :: calculate_rho0_atom, init_rho0
72 :
73 : CONTAINS
74 :
75 : ! **************************************************************************************************
76 : !> \brief ...
77 : !> \param mp_gau ...
78 : !> \param basis_1c ...
79 : !> \param harmonics ...
80 : !> \param nchannels ...
81 : !> \param nsotot ...
82 : ! **************************************************************************************************
83 3396 : SUBROUTINE calculate_mpole_gau(mp_gau, basis_1c, harmonics, nchannels, nsotot)
84 :
85 : TYPE(mpole_gau_overlap) :: mp_gau
86 : TYPE(gto_basis_set_type), POINTER :: basis_1c
87 : TYPE(harmonics_atom_type), POINTER :: harmonics
88 : INTEGER, INTENT(IN) :: nchannels, nsotot
89 :
90 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_mpole_gau'
91 :
92 : INTEGER :: handle, icg, ig1, ig2, ipgf1, ipgf2, iset1, iset2, iso, iso1, iso2, l, l1, l2, &
93 : llmax, m1, m2, max_iso_not0_local, max_s_harm, maxl, maxso, n1, n2, nset
94 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list
95 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list
96 3396 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf
97 : REAL(KIND=dp) :: zet1, zet2
98 3396 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zet
99 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: my_CG
100 :
101 3396 : CALL timeset(routineN, handle)
102 :
103 3396 : NULLIFY (lmax, lmin, npgf, my_CG, zet)
104 :
105 3396 : CALL reallocate(mp_gau%Qlm_gg, 1, nsotot, 1, nsotot, 1, nchannels)
106 :
107 : CALL get_gto_basis_set(gto_basis_set=basis_1c, &
108 : lmax=lmax, lmin=lmin, maxso=maxso, &
109 3396 : npgf=npgf, nset=nset, zet=zet, maxl=maxl)
110 :
111 3396 : max_s_harm = harmonics%max_s_harm
112 3396 : llmax = harmonics%llmax
113 :
114 20376 : ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm))
115 :
116 3396 : my_CG => harmonics%my_CG
117 :
118 3396 : m1 = 0
119 12774 : DO iset1 = 1, nset
120 : m2 = 0
121 43008 : DO iset2 = 1, nset
122 :
123 : CALL get_none0_cg_list(my_CG, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
124 33630 : max_s_harm, llmax, cg_list, cg_n_list, max_iso_not0_local)
125 :
126 33630 : n1 = nsoset(lmax(iset1))
127 114890 : DO ipgf1 = 1, npgf(iset1)
128 81260 : zet1 = zet(ipgf1, iset1)
129 :
130 81260 : n2 = nsoset(lmax(iset2))
131 334062 : DO ipgf2 = 1, npgf(iset2)
132 219172 : zet2 = zet(ipgf2, iset2)
133 :
134 1587224 : DO iso = 1, MIN(nchannels, max_iso_not0_local)
135 1286792 : l = indso(1, iso)
136 3703322 : DO icg = 1, cg_n_list(iso)
137 2197358 : iso1 = cg_list(1, icg, iso)
138 2197358 : iso2 = cg_list(2, icg, iso)
139 :
140 2197358 : l1 = indso(1, iso1)
141 2197358 : l2 = indso(1, iso2)
142 2197358 : ig1 = iso1 + n1*(ipgf1 - 1) + m1
143 2197358 : ig2 = iso2 + n2*(ipgf2 - 1) + m2
144 :
145 : mp_gau%Qlm_gg(ig1, ig2, iso) = fourpi/(2._dp*l + 1._dp)* &
146 3484150 : my_CG(iso1, iso2, iso)*gaussint_sph(zet1 + zet2, l + l1 + l2)
147 : END DO ! icg
148 : END DO ! iso
149 :
150 : END DO ! ipgf2
151 : END DO ! ipgf1
152 76638 : m2 = m2 + maxso
153 : END DO ! iset2
154 12774 : m1 = m1 + maxso
155 : END DO ! iset1
156 :
157 3396 : DEALLOCATE (cg_list, cg_n_list)
158 :
159 3396 : CALL timestop(handle)
160 6792 : END SUBROUTINE calculate_mpole_gau
161 :
162 : ! **************************************************************************************************
163 : !> \brief ...
164 : !> \param gapw_control ...
165 : !> \param rho_atom_set ...
166 : !> \param rho0_atom_set ...
167 : !> \param rho0_mp ...
168 : !> \param a_list ...
169 : !> \param natom ...
170 : !> \param ikind ...
171 : !> \param qs_kind ...
172 : !> \param rho0_h_tot ...
173 : ! **************************************************************************************************
174 32746 : SUBROUTINE calculate_rho0_atom(gapw_control, rho_atom_set, rho0_atom_set, &
175 32746 : rho0_mp, a_list, natom, ikind, qs_kind, rho0_h_tot)
176 :
177 : TYPE(gapw_control_type), POINTER :: gapw_control
178 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
179 : TYPE(rho0_atom_type), DIMENSION(:), POINTER :: rho0_atom_set
180 : TYPE(rho0_mpole_type), POINTER :: rho0_mp
181 : INTEGER, DIMENSION(:), INTENT(IN) :: a_list
182 : INTEGER, INTENT(IN) :: natom, ikind
183 : TYPE(qs_kind_type), INTENT(IN) :: qs_kind
184 : REAL(KIND=dp), INTENT(INOUT) :: rho0_h_tot
185 :
186 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_rho0_atom'
187 :
188 : INTEGER :: handle, iat, iatom, ic, ico, ir, is, &
189 : iso, ispin, l, lmax0, lshell, lx, ly, &
190 : lz, nr, nsotot, nspins
191 : LOGICAL :: paw_atom
192 : REAL(KIND=dp) :: sum1
193 32746 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: cpc_ah, cpc_as
194 32746 : REAL(KIND=dp), DIMENSION(:), POINTER :: norm_g0l_h
195 32746 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: g0_h, vg0_h
196 : TYPE(grid_atom_type), POINTER :: g_atom
197 : TYPE(harmonics_atom_type), POINTER :: harmonics
198 : TYPE(mpole_gau_overlap), POINTER :: mpole_gau
199 : TYPE(mpole_rho_atom), POINTER :: mpole_rho
200 32746 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: cpc_h, cpc_s
201 : TYPE(rho_atom_type), POINTER :: rho_atom
202 :
203 32746 : CALL timeset(routineN, handle)
204 :
205 32746 : NULLIFY (mpole_gau)
206 32746 : NULLIFY (mpole_rho)
207 32746 : NULLIFY (g0_h, vg0_h, g_atom)
208 32746 : NULLIFY (norm_g0l_h, harmonics)
209 :
210 : CALL get_rho0_mpole(rho0_mpole=rho0_mp, ikind=ikind, &
211 : l0_ikind=lmax0, mp_gau_ikind=mpole_gau, &
212 : g0_h=g0_h, &
213 : vg0_h=vg0_h, &
214 32746 : norm_g0l_h=norm_g0l_h)
215 :
216 32746 : CALL get_qs_kind(qs_kind, harmonics=harmonics, paw_atom=paw_atom, grid_atom=g_atom)
217 :
218 32746 : nr = g_atom%nr
219 :
220 : ! Set density coefficient to zero before the calculation
221 86642 : DO iat = 1, natom
222 53896 : iatom = a_list(iat)
223 20985872 : rho0_atom_set(iatom)%rho0_rad_h%r_coef = 0.0_dp
224 432992 : rho0_mp%mp_rho(iatom)%Qlm_tot = 0.0_dp
225 53896 : rho0_mp%mp_rho(iatom)%Qlm_tot(1) = rho0_mp%mp_rho(iatom)%Qlm_z
226 161688 : rho0_mp%mp_rho(iatom)%Q0 = 0.0_dp
227 510364 : rho0_mp%mp_rho(iatom)%Qlm_car = 0.0_dp
228 : END DO
229 :
230 32746 : IF (.NOT. (.NOT. paw_atom .AND. gapw_control%nopaw_as_gpw)) THEN
231 80532 : DO iat = 1, natom
232 49332 : iatom = a_list(iat)
233 49332 : mpole_rho => rho0_mp%mp_rho(iatom)
234 49332 : rho_atom => rho_atom_set(iatom)
235 :
236 49332 : IF (paw_atom) THEN
237 49332 : NULLIFY (cpc_h, cpc_s)
238 49332 : CALL get_rho_atom(rho_atom=rho_atom, cpc_h=cpc_h, cpc_s=cpc_s)
239 49332 : nspins = SIZE(cpc_h)
240 49332 : nsotot = SIZE(mpole_gau%Qlm_gg, 1)
241 246660 : ALLOCATE (cpc_ah(nsotot, nsotot, nspins))
242 161313096 : cpc_ah = 0._dp
243 246660 : ALLOCATE (cpc_as(nsotot, nsotot, nspins))
244 161313096 : cpc_as = 0._dp
245 106752 : DO ispin = 1, nspins
246 57420 : CALL prj_scatter(cpc_h(ispin)%r_coef, cpc_ah(:, :, ispin), qs_kind)
247 106752 : CALL prj_scatter(cpc_s(ispin)%r_coef, cpc_as(:, :, ispin), qs_kind)
248 : END DO
249 : END IF
250 :
251 : ! Total charge (hard-soft) at atom
252 49332 : IF (paw_atom) THEN
253 106752 : DO ispin = 1, nspins
254 : mpole_rho%Q0(ispin) = (trace_r_AxB(mpole_gau%Qlm_gg(:, :, 1), nsotot, &
255 : cpc_ah(:, :, ispin), nsotot, nsotot, nsotot) &
256 : - trace_r_AxB(mpole_gau%Qlm_gg(:, :, 1), nsotot, &
257 106752 : cpc_as(:, :, ispin), nsotot, nsotot, nsotot))/SQRT(fourpi)
258 : END DO
259 : END IF
260 : ! Multipoles of local charge distribution
261 423864 : DO iso = 1, nsoset(lmax0)
262 374532 : l = indso(1, iso)
263 374532 : IF (paw_atom) THEN
264 374532 : mpole_rho%Qlm_h(iso) = 0.0_dp
265 374532 : mpole_rho%Qlm_s(iso) = 0.0_dp
266 :
267 808224 : DO ispin = 1, nspins
268 : mpole_rho%Qlm_h(iso) = mpole_rho%Qlm_h(iso) + &
269 : trace_r_AxB(mpole_gau%Qlm_gg(:, :, iso), nsotot, &
270 433692 : cpc_ah(:, :, ispin), nsotot, nsotot, nsotot)
271 : mpole_rho%Qlm_s(iso) = mpole_rho%Qlm_s(iso) + &
272 : trace_r_AxB(mpole_gau%Qlm_gg(:, :, iso), nsotot, &
273 808224 : cpc_as(:, :, ispin), nsotot, nsotot, nsotot)
274 : END DO ! ispin
275 :
276 : mpole_rho%Qlm_tot(iso) = mpole_rho%Qlm_tot(iso) + &
277 374532 : mpole_rho%Qlm_h(iso) - mpole_rho%Qlm_s(iso)
278 : END IF
279 :
280 : rho0_atom_set(iatom)%rho0_rad_h%r_coef(1:nr, iso) = &
281 41024372 : g0_h(1:nr, l)*mpole_rho%Qlm_tot(iso)
282 : rho0_atom_set(iatom)%vrho0_rad_h%r_coef(1:nr, iso) = &
283 41024372 : vg0_h(1:nr, l)*mpole_rho%Qlm_tot(iso)
284 :
285 374532 : sum1 = 0.0_dp
286 20699452 : DO ir = 1, nr
287 : sum1 = sum1 + g_atom%wr(ir)* &
288 20699452 : rho0_atom_set(iatom)%rho0_rad_h%r_coef(ir, iso)
289 : END DO
290 423864 : rho0_h_tot = rho0_h_tot + sum1*harmonics%slm_int(iso)
291 : END DO ! iso
292 82078 : IF (paw_atom) THEN
293 49332 : DEALLOCATE (cpc_ah, cpc_as)
294 : END IF
295 : END DO ! iat
296 : END IF
297 :
298 : ! Transform the coefficinets from spherical to Cartesian
299 32746 : IF (.NOT. paw_atom .AND. gapw_control%nopaw_as_gpw) THEN
300 6110 : DO iat = 1, natom
301 4564 : iatom = a_list(iat)
302 4564 : mpole_rho => rho0_mp%mp_rho(iatom)
303 :
304 10674 : DO lshell = 0, lmax0
305 13692 : DO ic = 1, nco(lshell)
306 4564 : ico = ic + ncoset(lshell - 1)
307 9128 : mpole_rho%Qlm_car(ico) = 0.0_dp
308 : END DO
309 : END DO
310 : END DO
311 : ELSE
312 80532 : DO iat = 1, natom
313 49332 : iatom = a_list(iat)
314 49332 : mpole_rho => rho0_mp%mp_rho(iatom)
315 210028 : DO lshell = 0, lmax0
316 597986 : DO ic = 1, nco(lshell)
317 419158 : ico = ic + ncoset(lshell - 1)
318 419158 : mpole_rho%Qlm_car(ico) = 0.0_dp
319 419158 : lx = indco(1, ico)
320 419158 : ly = indco(2, ico)
321 419158 : lz = indco(3, ico)
322 2255472 : DO is = 1, nso(lshell)
323 1706818 : iso = is + nsoset(lshell - 1)
324 : mpole_rho%Qlm_car(ico) = mpole_rho%Qlm_car(ico) + &
325 : norm_g0l_h(lshell)* &
326 : orbtramat(lshell)%slm(is, ic)* &
327 2125976 : mpole_rho%Qlm_tot(iso)
328 :
329 : END DO
330 : END DO
331 : END DO ! lshell
332 : END DO ! iat
333 : END IF
334 : !MI Get rid of full gapw
335 :
336 32746 : CALL timestop(handle)
337 :
338 65492 : END SUBROUTINE calculate_rho0_atom
339 :
340 : ! **************************************************************************************************
341 : !> \brief ...
342 : !> \param local_rho_set ...
343 : !> \param qs_env ...
344 : !> \param gapw_control ...
345 : !> \param zcore ...
346 : ! **************************************************************************************************
347 3564 : SUBROUTINE init_rho0(local_rho_set, qs_env, gapw_control, zcore)
348 :
349 : TYPE(local_rho_type), POINTER :: local_rho_set
350 : TYPE(qs_environment_type), POINTER :: qs_env
351 : TYPE(gapw_control_type), POINTER :: gapw_control
352 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: zcore
353 :
354 : CHARACTER(len=*), PARAMETER :: routineN = 'init_rho0'
355 :
356 : CHARACTER(LEN=default_string_length) :: unit_str
357 : INTEGER :: handle, iat, iatom, ikind, l, l_rho1_max, laddg, lmaxg, maxl, maxnset, maxso, &
358 : nat, natom, nchan_c, nchan_s, nkind, nr, nset, nsotot, output_unit
359 1782 : INTEGER, DIMENSION(:), POINTER :: atom_list
360 : LOGICAL :: paw_atom
361 : REAL(KIND=dp) :: alpha_core, eps_Vrho0, max_rpgf0_s, &
362 : radius, rc_min, rc_orb, &
363 : total_rho_core_rspace, zeff
364 1782 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
365 : TYPE(cp_logger_type), POINTER :: logger
366 : TYPE(grid_atom_type), POINTER :: grid_atom
367 : TYPE(gto_basis_set_type), POINTER :: basis_1c
368 : TYPE(harmonics_atom_type), POINTER :: harmonics
369 1782 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
370 1782 : TYPE(rho0_atom_type), DIMENSION(:), POINTER :: rho0_atom_set
371 : TYPE(rho0_mpole_type), POINTER :: rho0_mpole
372 1782 : TYPE(rhoz_type), DIMENSION(:), POINTER :: rhoz_set
373 : TYPE(section_vals_type), POINTER :: dft_section
374 :
375 1782 : CALL timeset(routineN, handle)
376 :
377 1782 : NULLIFY (logger)
378 1782 : logger => cp_get_default_logger()
379 :
380 1782 : NULLIFY (qs_kind_set)
381 1782 : NULLIFY (atomic_kind_set)
382 1782 : NULLIFY (harmonics)
383 1782 : NULLIFY (basis_1c)
384 1782 : NULLIFY (rho0_mpole)
385 1782 : NULLIFY (rho0_atom_set)
386 1782 : NULLIFY (rhoz_set)
387 :
388 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
389 1782 : atomic_kind_set=atomic_kind_set)
390 :
391 1782 : nkind = SIZE(atomic_kind_set)
392 1782 : eps_Vrho0 = gapw_control%eps_Vrho0
393 :
394 : ! Initialize rhoz total to zero
395 : ! in gapw rhoz is calculated on local the lebedev grids
396 1782 : total_rho_core_rspace = 0.0_dp
397 :
398 1782 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
399 :
400 : ! Initialize the multipole and the compensation charge type
401 1782 : CALL allocate_rho0_mpole(rho0_mpole)
402 1782 : CALL allocate_rho0_atom(rho0_atom_set, natom)
403 :
404 : ! Allocate the multipole set
405 1782 : CALL allocate_multipoles(rho0_mpole%mp_rho, natom, rho0_mpole%mp_gau, nkind)
406 :
407 : ! Allocate the core density on the radial grid for each kind: rhoz_set
408 1782 : CALL allocate_rhoz(rhoz_set, nkind)
409 :
410 : ! For each kind, determine the max l for the compensation charge density
411 1782 : lmaxg = gapw_control%lmax_rho0
412 1782 : laddg = gapw_control%ladd_rho0
413 :
414 1782 : CALL reallocate(rho0_mpole%lmax0_kind, 1, nkind)
415 :
416 1782 : rho0_mpole%lmax_0 = 0
417 1782 : rc_min = 100.0_dp
418 1782 : maxnset = 0
419 5520 : DO ikind = 1, nkind
420 3738 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
421 : CALL get_qs_kind(qs_kind_set(ikind), &
422 : ngrid_rad=nr, &
423 : grid_atom=grid_atom, &
424 : harmonics=harmonics, &
425 : paw_atom=paw_atom, &
426 : hard0_radius=rc_orb, &
427 : zeff=zeff, &
428 3738 : alpha_core_charge=alpha_core)
429 : CALL get_qs_kind(qs_kind_set(ikind), &
430 3738 : basis_set=basis_1c, basis_type="GAPW_1C")
431 :
432 : ! Set charge distribution of ionic cores to zero when computing the response-density
433 3738 : IF (PRESENT(zcore)) zeff = zcore
434 :
435 : CALL get_gto_basis_set(gto_basis_set=basis_1c, &
436 : maxl=maxl, &
437 3738 : maxso=maxso, nset=nset)
438 :
439 3738 : maxnset = MAX(maxnset, nset)
440 :
441 3738 : l_rho1_max = indso(1, harmonics%max_iso_not0)
442 3738 : IF (paw_atom) THEN
443 3396 : rho0_mpole%lmax0_kind(ikind) = MIN(2*maxl, l_rho1_max, maxl + laddg, lmaxg)
444 : ELSE
445 342 : rho0_mpole%lmax0_kind(ikind) = 0
446 : END IF
447 :
448 3738 : CALL set_qs_kind(qs_kind_set(ikind), lmax_rho0=rho0_mpole%lmax0_kind(ikind))
449 :
450 3738 : IF (gapw_control%lrho1_eq_lrho0) harmonics%max_iso_not0 = &
451 0 : nsoset(rho0_mpole%lmax0_kind(ikind))
452 :
453 3738 : rho0_mpole%lmax_0 = MAX(rho0_mpole%lmax_0, rho0_mpole%lmax0_kind(ikind))
454 3738 : rc_min = MIN(rc_min, rc_orb)
455 :
456 3738 : nchan_s = nsoset(rho0_mpole%lmax0_kind(ikind))
457 3738 : nchan_c = ncoset(rho0_mpole%lmax0_kind(ikind))
458 3738 : nsotot = maxso*nset
459 :
460 9856 : DO iat = 1, nat
461 6118 : iatom = atom_list(iat)
462 : ! Allocate the multipole for rho1_h rho1_s and rho_z
463 6118 : CALL initialize_mpole_rho(rho0_mpole%mp_rho(iatom), nchan_s, nchan_c, zeff)
464 : ! Allocate the radial part of rho0_h and rho0_s
465 : ! This is calculated on the radial grid centered at the atomic position
466 9856 : CALL allocate_rho0_atom_rad(rho0_atom_set(iatom), nr, nchan_s)
467 : END DO
468 :
469 3738 : IF (paw_atom) THEN
470 : ! Calculate multipoles given by the product of 2 primitives Qlm_gg
471 : CALL calculate_mpole_gau(rho0_mpole%mp_gau(ikind), &
472 3396 : basis_1c, harmonics, nchan_s, nsotot)
473 : END IF
474 :
475 : ! Calculate the core density rhoz
476 : ! exp(-alpha_c**2 r**2)Z(alpha_c**2/pi)**(3/2)
477 : ! on the logarithmic radial grid
478 : ! WARNING: alpha_core_charge = alpha_c**2
479 : CALL calculate_rhoz(rhoz_set(ikind), grid_atom, alpha_core, zeff, &
480 12996 : nat, total_rho_core_rspace, harmonics)
481 : END DO ! ikind
482 1782 : total_rho_core_rspace = -total_rho_core_rspace
483 :
484 1782 : IF (gapw_control%alpha0_hard_from_input) THEN
485 : ! The exponent for the compensation charge rho0_hard is read from input
486 104 : rho0_mpole%zet0_h = gapw_control%alpha0_hard
487 : ELSE
488 : ! Calculate the exponent for the compensation charge rho0_hard
489 1678 : rho0_mpole%zet0_h = 0.1_dp
490 154834 : DO
491 156512 : radius = exp_radius(rho0_mpole%lmax_0, rho0_mpole%zet0_h, eps_Vrho0, 1.0_dp)
492 156512 : IF (radius <= rc_min) EXIT
493 154834 : rho0_mpole%zet0_h = rho0_mpole%zet0_h + 0.1_dp
494 : END DO
495 :
496 : END IF
497 :
498 : ! Allocate and calculate the normalization factors for g0_lm_h and g0_lm_s
499 1782 : CALL reallocate(rho0_mpole%norm_g0l_h, 0, rho0_mpole%lmax_0)
500 7140 : DO l = 0, rho0_mpole%lmax_0
501 : rho0_mpole%norm_g0l_h(l) = (2._dp*l + 1._dp)/ &
502 7140 : (fourpi*gaussint_sph(rho0_mpole%zet0_h, 2*l))
503 : END DO
504 :
505 : ! Allocate and Initialize the g0 gaussians used to build the compensation density
506 : ! and calculate the interaction radii
507 1782 : max_rpgf0_s = 0.0_dp
508 5520 : DO ikind = 1, nkind
509 3738 : CALL get_qs_kind(qs_kind_set(ikind), grid_atom=grid_atom)
510 3738 : CALL calculate_g0(rho0_mpole, grid_atom, ikind)
511 5520 : CALL interaction_radii_g0(rho0_mpole, ikind, eps_Vrho0, max_rpgf0_s)
512 : END DO
513 1782 : rho0_mpole%max_rpgf0_s = max_rpgf0_s
514 :
515 1782 : CALL set_local_rho(local_rho_set, rho0_atom_set=rho0_atom_set, rho0_mpole=rho0_mpole, rhoz_set=rhoz_set)
516 1782 : local_rho_set%rhoz_tot = total_rho_core_rspace
517 :
518 1782 : dft_section => section_vals_get_subs_vals(qs_env%input, "DFT")
519 : output_unit = cp_print_key_unit_nr(logger, dft_section, "PRINT%GAPW%RHO0_INFORMATION", &
520 1782 : extension=".Log")
521 1782 : CALL section_vals_val_get(dft_section, "PRINT%GAPW%RHO0_INFORMATION%UNIT", c_val=unit_str)
522 1782 : IF (output_unit > 0) THEN
523 1 : CALL write_rho0_info(rho0_mpole, unit_str, output_unit)
524 : END IF
525 : CALL cp_print_key_finished_output(output_unit, logger, dft_section, &
526 1782 : "PRINT%GAPW%RHO0_INFORMATION")
527 :
528 1782 : CALL timestop(handle)
529 :
530 1782 : END SUBROUTINE init_rho0
531 :
532 : ! **************************************************************************************************
533 : !> \brief ...
534 : !> \param rho0_mpole ...
535 : !> \param ik ...
536 : !> \param eps_Vrho0 ...
537 : !> \param max_rpgf0_s ...
538 : ! **************************************************************************************************
539 3738 : SUBROUTINE interaction_radii_g0(rho0_mpole, ik, eps_Vrho0, max_rpgf0_s)
540 :
541 : TYPE(rho0_mpole_type), POINTER :: rho0_mpole
542 : INTEGER, INTENT(IN) :: ik
543 : REAL(KIND=dp), INTENT(IN) :: eps_Vrho0
544 : REAL(KIND=dp), INTENT(INOUT) :: max_rpgf0_s
545 :
546 : INTEGER :: l, lmax
547 : REAL(KIND=dp) :: r_h, z0_h
548 3738 : REAL(KIND=dp), DIMENSION(:), POINTER :: ng0_h
549 :
550 : CALL get_rho0_mpole(rho0_mpole, ikind=ik, l0_ikind=lmax, &
551 3738 : zet0_h=z0_h, norm_g0l_h=ng0_h)
552 3738 : r_h = 0.0_dp
553 13976 : DO l = 0, lmax
554 13976 : r_h = MAX(r_h, exp_radius(l, z0_h, eps_Vrho0, ng0_h(l), rlow=r_h))
555 : END DO
556 :
557 3738 : rho0_mpole%mp_gau(ik)%rpgf0_h = r_h
558 3738 : rho0_mpole%mp_gau(ik)%rpgf0_s = r_h
559 3738 : max_rpgf0_s = MAX(max_rpgf0_s, r_h)
560 :
561 3738 : END SUBROUTINE interaction_radii_g0
562 :
563 : END MODULE qs_rho0_methods
|
