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 : !> \brief Calculate the interaction radii for the operator matrix calculation.
10 : !> \par History
11 : !> Joost VandeVondele : added exp_radius_very_extended
12 : !> 24.09.2002 overloaded init_interaction_radii for KG use (gt)
13 : !> 27.02.2004 integrated KG into QS version of init_interaction_radii (jgh)
14 : !> 07.03.2006 new routine to extend kind interaction radius for semi-empiric (jgh)
15 : !> \author MK (12.07.2000)
16 : ! **************************************************************************************************
17 : MODULE qs_interactions
18 :
19 : USE ao_util, ONLY: exp_radius
20 : USE atomic_kind_types, ONLY: atomic_kind_type,&
21 : get_atomic_kind
22 : USE basis_set_types, ONLY: get_gto_basis_set,&
23 : gto_basis_set_type,&
24 : set_gto_basis_set
25 : USE cp_control_types, ONLY: qs_control_type,&
26 : semi_empirical_control_type
27 : USE cp_log_handling, ONLY: cp_get_default_logger,&
28 : cp_logger_type
29 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
30 : cp_print_key_unit_nr
31 : USE cp_units, ONLY: cp_unit_from_cp2k
32 : USE external_potential_types, ONLY: all_potential_type,&
33 : get_potential,&
34 : gth_potential_type,&
35 : set_potential,&
36 : sgp_potential_type
37 : USE input_section_types, ONLY: section_vals_type,&
38 : section_vals_val_get
39 : USE kinds, ONLY: default_string_length,&
40 : dp
41 : USE paw_proj_set_types, ONLY: get_paw_proj_set,&
42 : paw_proj_set_type,&
43 : set_paw_proj_set
44 : USE qs_kind_types, ONLY: get_qs_kind,&
45 : qs_kind_type
46 : USE string_utilities, ONLY: uppercase
47 : #include "./base/base_uses.f90"
48 :
49 : IMPLICIT NONE
50 :
51 : PRIVATE
52 :
53 : ! *** Global parameters (only in this module)
54 :
55 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_interactions'
56 :
57 : ! *** Public subroutines ***
58 :
59 : PUBLIC :: init_interaction_radii, init_se_nlradius, init_interaction_radii_orb_basis
60 :
61 : PUBLIC :: write_paw_radii, write_ppnl_radii, write_ppl_radii, write_core_charge_radii, &
62 : write_pgf_orb_radii
63 :
64 : CONTAINS
65 :
66 : ! **************************************************************************************************
67 : !> \brief Initialize all the atomic kind radii for a given threshold value.
68 : !> \param qs_control ...
69 : !> \param qs_kind_set ...
70 : !> \date 24.09.2002
71 : !> \author GT
72 : !> \version 1.0
73 : ! **************************************************************************************************
74 8556 : SUBROUTINE init_interaction_radii(qs_control, qs_kind_set)
75 :
76 : TYPE(qs_control_type), INTENT(IN) :: qs_control
77 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
78 :
79 : CHARACTER(len=*), PARAMETER :: routineN = 'init_interaction_radii'
80 :
81 : INTEGER :: handle, i, iexp_ppl, ikind, ip, iprj_ppnl, j, l, lppl, lppnl, lppsl, lprj, &
82 : lprj_ppnl, maxl, n_local, nexp_lpot, nexp_lsd, nexp_ppl, nkind, nloc
83 : INTEGER, DIMENSION(0:10) :: npot
84 : INTEGER, DIMENSION(1:10) :: nrloc
85 : INTEGER, DIMENSION(1:15, 0:10) :: nrpot
86 8556 : INTEGER, DIMENSION(:), POINTER :: nct_lpot, nct_lsd, nprj, nprj_ppnl
87 : LOGICAL :: ecp_local, ecp_semi_local, llsd, lpot, &
88 : paw_atom
89 : LOGICAL, DIMENSION(0:5) :: is_nonlocal
90 : REAL(KIND=dp) :: alpha_core_charge, alpha_ppl, ccore_charge, cerf_ppl, core_charge_radius, &
91 : cval, ppl_radius, ppnl_radius, rcprj, zeta
92 : REAL(KIND=dp), DIMENSION(1:10) :: aloc, bloc
93 : REAL(KIND=dp), DIMENSION(1:15, 0:10) :: apot, bpot
94 8556 : REAL(KIND=dp), DIMENSION(:), POINTER :: a_local, a_nonlocal, alpha_lpot, &
95 8556 : alpha_lsd, alpha_ppnl, c_local, &
96 8556 : cexp_ppl
97 8556 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cprj_ppnl, cval_lpot, cval_lsd, rzetprj, &
98 8556 : zet
99 8556 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: c_nonlocal
100 : TYPE(all_potential_type), POINTER :: all_potential
101 : TYPE(gth_potential_type), POINTER :: gth_potential
102 : TYPE(gto_basis_set_type), POINTER :: aux_basis_set, aux_fit_basis_set, aux_gw_basis, &
103 : aux_opt_basis_set, gapw_1c_basis, harris_basis, lri_basis, mao_basis, min_basis_set, &
104 : orb_basis_set, p_lri_basis, ri_aux_basis_set, ri_basis, ri_xas_basis, soft_basis, &
105 : tda_k_basis
106 : TYPE(paw_proj_set_type), POINTER :: paw_proj_set
107 : TYPE(sgp_potential_type), POINTER :: sgp_potential
108 :
109 8556 : CALL timeset(routineN, handle)
110 :
111 8556 : NULLIFY (all_potential, gth_potential, sgp_potential)
112 8556 : NULLIFY (aux_basis_set, aux_fit_basis_set, aux_gw_basis, tda_k_basis, &
113 8556 : harris_basis, lri_basis, mao_basis, orb_basis_set, p_lri_basis, ri_aux_basis_set, &
114 8556 : ri_basis, ri_xas_basis, soft_basis, gapw_1c_basis, aux_opt_basis_set, min_basis_set)
115 8556 : NULLIFY (nprj_ppnl, nprj)
116 8556 : NULLIFY (alpha_ppnl, cexp_ppl, cprj_ppnl, zet)
117 :
118 8556 : nkind = SIZE(qs_kind_set)
119 8556 : nexp_ppl = 0
120 :
121 24722 : DO ikind = 1, nkind
122 :
123 16166 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type="ORB")
124 16166 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_basis_set, basis_type="AUX")
125 16166 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=min_basis_set, basis_type="MIN")
126 16166 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_fit_basis_set, basis_type="AUX_FIT")
127 16166 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_opt_basis_set, basis_type="AUX_OPT")
128 16166 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=lri_basis, basis_type="LRI_AUX")
129 16166 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=p_lri_basis, basis_type="P_LRI_AUX")
130 16166 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=ri_basis, basis_type="RI_HXC")
131 16166 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=ri_aux_basis_set, basis_type="RI_AUX")
132 16166 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=mao_basis, basis_type="MAO")
133 16166 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=harris_basis, basis_type="HARRIS")
134 16166 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_gw_basis, basis_type="AUX_GW")
135 16166 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=ri_xas_basis, basis_type="RI_XAS")
136 16166 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=soft_basis, basis_type="ORB_SOFT")
137 16166 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=gapw_1c_basis, basis_type="GAPW_1C")
138 16166 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=tda_k_basis, basis_type="TDA_HFX")
139 : CALL get_qs_kind(qs_kind_set(ikind), &
140 : paw_proj_set=paw_proj_set, &
141 : paw_atom=paw_atom, &
142 : all_potential=all_potential, &
143 : gth_potential=gth_potential, &
144 16166 : sgp_potential=sgp_potential)
145 :
146 : ! Calculate the orbital basis function radii ***
147 : ! For ALL electron this has to come before the calculation of the
148 : ! radii used to calculate the 3c lists.
149 : ! Indeed, there the radii of the GTO primitives are needed
150 16166 : IF (ASSOCIATED(orb_basis_set)) THEN
151 : CALL init_interaction_radii_orb_basis(orb_basis_set, qs_control%eps_pgf_orb, &
152 15684 : qs_control%eps_kg_orb)
153 : END IF
154 :
155 16166 : IF (ASSOCIATED(all_potential)) THEN
156 :
157 : CALL get_potential(potential=all_potential, &
158 : alpha_core_charge=alpha_core_charge, &
159 4370 : ccore_charge=ccore_charge)
160 :
161 : ! Calculate the radius of the core charge distribution
162 :
163 : core_charge_radius = exp_radius(0, alpha_core_charge, &
164 : qs_control%eps_core_charge, &
165 4370 : ccore_charge)
166 :
167 : CALL set_potential(potential=all_potential, &
168 4370 : core_charge_radius=core_charge_radius)
169 :
170 11796 : ELSE IF (ASSOCIATED(gth_potential)) THEN
171 :
172 : CALL get_potential(potential=gth_potential, &
173 : alpha_core_charge=alpha_core_charge, &
174 : ccore_charge=ccore_charge, &
175 : alpha_ppl=alpha_ppl, &
176 : cerf_ppl=cerf_ppl, &
177 : nexp_ppl=nexp_ppl, &
178 : cexp_ppl=cexp_ppl, &
179 : lpot_present=lpot, &
180 : lsd_present=llsd, &
181 : lppnl=lppnl, &
182 : alpha_ppnl=alpha_ppnl, &
183 : nprj_ppnl=nprj_ppnl, &
184 11604 : cprj_ppnl=cprj_ppnl)
185 :
186 : ! Calculate the radius of the core charge distribution
187 : core_charge_radius = exp_radius(0, alpha_core_charge, &
188 : qs_control%eps_core_charge, &
189 11604 : ccore_charge)
190 :
191 : ! Calculate the radii of the local part
192 : ! of the Goedecker pseudopotential (GTH)
193 11604 : ppl_radius = exp_radius(0, alpha_ppl, qs_control%eps_ppl, cerf_ppl)
194 :
195 34372 : DO iexp_ppl = 1, nexp_ppl
196 22768 : lppl = 2*(iexp_ppl - 1)
197 : ppl_radius = MAX(ppl_radius, &
198 : exp_radius(lppl, alpha_ppl, &
199 : qs_control%eps_ppl, &
200 : cexp_ppl(iexp_ppl), &
201 34372 : rlow=ppl_radius))
202 : END DO
203 11604 : IF (lpot) THEN
204 : ! extended local potential
205 : CALL get_potential(potential=gth_potential, &
206 : nexp_lpot=nexp_lpot, &
207 : alpha_lpot=alpha_lpot, &
208 : nct_lpot=nct_lpot, &
209 8 : cval_lpot=cval_lpot)
210 20 : DO j = 1, nexp_lpot
211 46 : DO i = 1, nct_lpot(j)
212 26 : lppl = 2*(i - 1)
213 : ppl_radius = MAX(ppl_radius, &
214 : exp_radius(lppl, alpha_lpot(j), qs_control%eps_ppl, &
215 38 : cval_lpot(i, j), rlow=ppl_radius))
216 : END DO
217 : END DO
218 : END IF
219 11604 : IF (llsd) THEN
220 : ! lsd local potential
221 : CALL get_potential(potential=gth_potential, &
222 : nexp_lsd=nexp_lsd, &
223 : alpha_lsd=alpha_lsd, &
224 : nct_lsd=nct_lsd, &
225 0 : cval_lsd=cval_lsd)
226 0 : DO j = 1, nexp_lsd
227 0 : DO i = 1, nct_lsd(j)
228 0 : lppl = 2*(i - 1)
229 : ppl_radius = MAX(ppl_radius, &
230 : exp_radius(lppl, alpha_lsd(j), qs_control%eps_ppl, &
231 0 : cval_lsd(i, j), rlow=ppl_radius))
232 : END DO
233 : END DO
234 : END IF
235 :
236 : ! Calculate the radii of the non-local part
237 : ! of the Goedecker pseudopotential (GTH)
238 11604 : ppnl_radius = 0.0_dp
239 22120 : DO l = 0, lppnl
240 29494 : DO iprj_ppnl = 1, nprj_ppnl(l)
241 7374 : lprj_ppnl = l + 2*(iprj_ppnl - 1)
242 : ppnl_radius = MAX(ppnl_radius, &
243 : exp_radius(lprj_ppnl, alpha_ppnl(l), &
244 : qs_control%eps_pgf_orb, &
245 : cprj_ppnl(iprj_ppnl, l), &
246 17890 : rlow=ppnl_radius))
247 : END DO
248 : END DO
249 : CALL set_potential(potential=gth_potential, &
250 : core_charge_radius=core_charge_radius, &
251 : ppl_radius=ppl_radius, &
252 11604 : ppnl_radius=ppnl_radius)
253 :
254 192 : ELSE IF (ASSOCIATED(sgp_potential)) THEN
255 :
256 : ! Calculate the radius of the core charge distribution
257 : CALL get_potential(potential=sgp_potential, &
258 : alpha_core_charge=alpha_core_charge, &
259 28 : ccore_charge=ccore_charge)
260 : core_charge_radius = exp_radius(0, alpha_core_charge, &
261 : qs_control%eps_core_charge, &
262 28 : ccore_charge)
263 : ! Calculate the radii of the local part
264 28 : ppl_radius = core_charge_radius
265 28 : CALL get_potential(potential=sgp_potential, ecp_local=ecp_local)
266 28 : IF (ecp_local) THEN
267 16 : CALL get_potential(potential=sgp_potential, nloc=nloc, nrloc=nrloc, aloc=aloc, bloc=bloc)
268 68 : DO i = 1, nloc
269 52 : lppl = MAX(0, nrloc(i) - 2)
270 : ppl_radius = MAX(ppl_radius, &
271 68 : exp_radius(lppl, bloc(i), qs_control%eps_ppl, aloc(i), rlow=ppl_radius))
272 : END DO
273 : ELSE
274 12 : CALL get_potential(potential=sgp_potential, n_local=n_local, a_local=a_local, c_local=c_local)
275 156 : DO i = 1, n_local
276 : ppl_radius = MAX(ppl_radius, &
277 156 : exp_radius(0, a_local(i), qs_control%eps_ppl, c_local(i), rlow=ppl_radius))
278 : END DO
279 : END IF
280 : ! Calculate the radii of the semi-local part
281 28 : CALL get_potential(potential=sgp_potential, ecp_semi_local=ecp_semi_local)
282 28 : IF (ecp_semi_local) THEN
283 : CALL get_potential(potential=sgp_potential, sl_lmax=lppsl, npot=npot, nrpot=nrpot, &
284 16 : apot=apot, bpot=bpot)
285 68 : DO l = 0, lppsl
286 272 : DO i = 1, npot(l)
287 204 : lppl = MAX(0, nrpot(i, l) - 2)
288 : ppl_radius = MAX(ppl_radius, &
289 : exp_radius(lppl, bpot(i, l), qs_control%eps_ppl, apot(i, l), &
290 256 : rlow=ppl_radius))
291 : END DO
292 : END DO
293 : END IF
294 : ! Calculate the radii of the non-local part
295 28 : ppnl_radius = 0.0_dp
296 28 : CALL get_potential(potential=sgp_potential, lmax=lppnl, n_nonlocal=nloc)
297 : CALL get_potential(potential=sgp_potential, is_nonlocal=is_nonlocal, &
298 28 : a_nonlocal=a_nonlocal, c_nonlocal=c_nonlocal)
299 34 : DO l = 0, lppnl
300 34 : IF (is_nonlocal(l)) THEN
301 54 : DO i = 1, nloc
302 480 : cval = MAXVAL(ABS(c_nonlocal(i, :, l)))
303 : ppnl_radius = MAX(ppnl_radius, &
304 54 : exp_radius(l, a_nonlocal(i), qs_control%eps_pgf_orb, cval, rlow=ppnl_radius))
305 : END DO
306 : END IF
307 : END DO
308 : CALL set_potential(potential=sgp_potential, &
309 : core_charge_radius=core_charge_radius, &
310 : ppl_radius=ppl_radius, &
311 28 : ppnl_radius=ppnl_radius)
312 : END IF
313 :
314 : ! Calculate the aux fit orbital basis function radii
315 16166 : IF (ASSOCIATED(aux_fit_basis_set)) THEN
316 : CALL init_interaction_radii_orb_basis(aux_fit_basis_set, qs_control%eps_pgf_orb, &
317 1224 : qs_control%eps_kg_orb)
318 : END IF
319 :
320 : ! Calculate the aux opt orbital basis function radii
321 16166 : IF (ASSOCIATED(aux_opt_basis_set)) THEN
322 420 : CALL init_interaction_radii_orb_basis(aux_opt_basis_set, qs_control%eps_pgf_orb)
323 : END IF
324 :
325 : ! Calculate the minimal basis function radii
326 16166 : IF (ASSOCIATED(min_basis_set)) THEN
327 4 : CALL init_interaction_radii_orb_basis(min_basis_set, qs_control%eps_pgf_orb)
328 : END IF
329 :
330 : ! Calculate the aux orbital basis function radii
331 16166 : IF (ASSOCIATED(aux_basis_set)) THEN
332 0 : CALL init_interaction_radii_orb_basis(aux_basis_set, qs_control%eps_pgf_orb)
333 : END IF
334 :
335 : ! Calculate the RI aux orbital basis function radii
336 16166 : IF (ASSOCIATED(RI_aux_basis_set)) THEN
337 3926 : CALL init_interaction_radii_orb_basis(RI_aux_basis_set, qs_control%eps_pgf_orb)
338 : END IF
339 :
340 : ! Calculate the MAO orbital basis function radii
341 16166 : IF (ASSOCIATED(mao_basis)) THEN
342 20 : CALL init_interaction_radii_orb_basis(mao_basis, qs_control%eps_pgf_orb)
343 : END IF
344 :
345 : ! Calculate the HARRIS orbital basis function radii
346 16166 : IF (ASSOCIATED(harris_basis)) THEN
347 132 : CALL init_interaction_radii_orb_basis(harris_basis, qs_control%eps_pgf_orb)
348 : END IF
349 :
350 : ! Calculate the AUX_GW orbital basis function radii
351 16166 : IF (ASSOCIATED(aux_gw_basis)) THEN
352 36 : CALL init_interaction_radii_orb_basis(aux_gw_basis, qs_control%eps_pgf_orb)
353 : END IF
354 :
355 : ! Calculate the soft orbital basis function radii
356 16166 : IF (ASSOCIATED(soft_basis)) THEN
357 1952 : CALL init_interaction_radii_orb_basis(soft_basis, qs_control%eps_pgf_orb)
358 : END IF
359 :
360 : ! Calculate the GAPW 1 center basis function radii
361 16166 : IF (ASSOCIATED(gapw_1c_basis)) THEN
362 1952 : CALL init_interaction_radii_orb_basis(gapw_1c_basis, qs_control%eps_pgf_orb)
363 : END IF
364 :
365 : ! Calculate the HFX SR basis for TDA
366 16166 : IF (ASSOCIATED(tda_k_basis)) THEN
367 0 : CALL init_interaction_radii_orb_basis(tda_k_basis, qs_control%eps_pgf_orb)
368 : END IF
369 :
370 : ! Calculate the lri basis function radii
371 16166 : IF (ASSOCIATED(lri_basis)) THEN
372 86 : CALL init_interaction_radii_orb_basis(lri_basis, qs_control%eps_pgf_orb)
373 : END IF
374 16166 : IF (ASSOCIATED(ri_basis)) THEN
375 0 : CALL init_interaction_radii_orb_basis(ri_basis, qs_control%eps_pgf_orb)
376 : END IF
377 16166 : IF (ASSOCIATED(ri_xas_basis)) THEN
378 78 : CALL init_interaction_radii_orb_basis(ri_xas_basis, qs_control%eps_pgf_orb)
379 : END IF
380 16166 : IF (ASSOCIATED(p_lri_basis)) THEN
381 8 : CALL init_interaction_radii_orb_basis(p_lri_basis, qs_control%eps_pgf_orb)
382 : END IF
383 :
384 : ! Calculate the paw projector basis function radii
385 40888 : IF (ASSOCIATED(paw_proj_set)) THEN
386 : CALL get_paw_proj_set(paw_proj_set=paw_proj_set, &
387 : nprj=nprj, &
388 : maxl=maxl, &
389 : zetprj=zet, &
390 1642 : rzetprj=rzetprj)
391 1642 : rcprj = 0.0_dp
392 31526 : rzetprj = 0.0_dp
393 5554 : DO lprj = 0, maxl
394 22132 : DO ip = 1, nprj(lprj)
395 16578 : zeta = zet(ip, lprj)
396 : rzetprj(ip, lprj) = MAX(rzetprj(ip, lprj), &
397 : exp_radius(lprj, zeta, qs_control%eps_pgf_orb, &
398 16578 : 1.0_dp, rlow=rzetprj(ip, lprj)))
399 20490 : rcprj = MAX(rcprj, rzetprj(ip, lprj))
400 : END DO
401 : END DO
402 : CALL set_paw_proj_set(paw_proj_set=paw_proj_set, &
403 : rzetprj=rzetprj, &
404 1642 : rcprj=rcprj)
405 :
406 : END IF
407 : END DO ! ikind
408 :
409 8556 : CALL timestop(handle)
410 :
411 8556 : END SUBROUTINE init_interaction_radii
412 :
413 : ! **************************************************************************************************
414 : !> \brief ...
415 : !> \param orb_basis_set ...
416 : !> \param eps_pgf_orb ...
417 : !> \param eps_pgf_short ...
418 : ! **************************************************************************************************
419 29752 : SUBROUTINE init_interaction_radii_orb_basis(orb_basis_set, eps_pgf_orb, eps_pgf_short)
420 :
421 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
422 : REAL(KIND=dp), INTENT(IN) :: eps_pgf_orb
423 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: eps_pgf_short
424 :
425 : INTEGER :: ipgf, iset, ishell, l, nset
426 29752 : INTEGER, DIMENSION(:), POINTER :: npgf, nshell
427 29752 : INTEGER, DIMENSION(:, :), POINTER :: lshell
428 : REAL(KIND=dp) :: eps_short, gcca, kind_radius, &
429 : short_radius, zeta
430 29752 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius
431 29752 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pgf_radius, zet
432 29752 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: gcc
433 :
434 59504 : IF (ASSOCIATED(orb_basis_set)) THEN
435 :
436 29752 : IF (PRESENT(eps_pgf_short)) THEN
437 16908 : eps_short = eps_pgf_short
438 : ELSE
439 12844 : eps_short = eps_pgf_orb
440 : END IF
441 :
442 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
443 : nset=nset, &
444 : nshell=nshell, &
445 : npgf=npgf, &
446 : l=lshell, &
447 : zet=zet, &
448 : gcc=gcc, &
449 : pgf_radius=pgf_radius, &
450 29752 : set_radius=set_radius)
451 :
452 29752 : kind_radius = 0.0_dp
453 29752 : short_radius = 0.0_dp
454 :
455 141022 : DO iset = 1, nset
456 111270 : set_radius(iset) = 0.0_dp
457 327830 : DO ipgf = 1, npgf(iset)
458 216560 : pgf_radius(ipgf, iset) = 0.0_dp
459 731071 : DO ishell = 1, nshell(iset)
460 514511 : l = lshell(ishell, iset)
461 514511 : gcca = gcc(ipgf, ishell, iset)
462 514511 : zeta = zet(ipgf, iset)
463 : pgf_radius(ipgf, iset) = MAX(pgf_radius(ipgf, iset), &
464 : exp_radius(l, zeta, eps_pgf_orb, gcca, &
465 514511 : rlow=pgf_radius(ipgf, iset)))
466 : short_radius = MAX(short_radius, &
467 : exp_radius(l, zeta, eps_short, gcca, &
468 731071 : rlow=short_radius))
469 : END DO
470 327830 : set_radius(iset) = MAX(set_radius(iset), pgf_radius(ipgf, iset))
471 : END DO
472 141022 : kind_radius = MAX(kind_radius, set_radius(iset))
473 : END DO
474 :
475 : CALL set_gto_basis_set(gto_basis_set=orb_basis_set, &
476 : pgf_radius=pgf_radius, &
477 : set_radius=set_radius, &
478 : kind_radius=kind_radius, &
479 29752 : short_kind_radius=short_radius)
480 :
481 : END IF
482 :
483 29752 : END SUBROUTINE init_interaction_radii_orb_basis
484 :
485 : ! **************************************************************************************************
486 : !> \brief ...
487 : !> \param se_control ...
488 : !> \param atomic_kind_set ...
489 : !> \param qs_kind_set ...
490 : !> \param subsys_section ...
491 : ! **************************************************************************************************
492 996 : SUBROUTINE init_se_nlradius(se_control, atomic_kind_set, qs_kind_set, &
493 : subsys_section)
494 :
495 : TYPE(semi_empirical_control_type), POINTER :: se_control
496 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
497 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
498 : TYPE(section_vals_type), POINTER :: subsys_section
499 :
500 : INTEGER :: ikind, nkind
501 : REAL(KIND=dp) :: kind_radius
502 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
503 : TYPE(qs_kind_type), POINTER :: qs_kind
504 :
505 996 : NULLIFY (orb_basis_set, qs_kind)
506 :
507 996 : nkind = SIZE(qs_kind_set)
508 :
509 3232 : DO ikind = 1, nkind
510 :
511 2236 : qs_kind => qs_kind_set(ikind)
512 :
513 2236 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis_set)
514 :
515 3232 : IF (ASSOCIATED(orb_basis_set)) THEN
516 :
517 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
518 2236 : kind_radius=kind_radius)
519 :
520 2236 : kind_radius = MAX(kind_radius, se_control%cutoff_exc)
521 :
522 : CALL set_gto_basis_set(gto_basis_set=orb_basis_set, &
523 2236 : kind_radius=kind_radius)
524 :
525 : END IF
526 :
527 : END DO
528 :
529 996 : CALL write_kind_radii(atomic_kind_set, qs_kind_set, subsys_section)
530 :
531 996 : END SUBROUTINE init_se_nlradius
532 :
533 : ! **************************************************************************************************
534 : !> \brief ...
535 : !> \param atomic_kind_set ...
536 : !> \param qs_kind_set ...
537 : !> \param subsys_section ...
538 : ! **************************************************************************************************
539 996 : SUBROUTINE write_kind_radii(atomic_kind_set, qs_kind_set, subsys_section)
540 :
541 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
542 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
543 : TYPE(section_vals_type), POINTER :: subsys_section
544 :
545 : CHARACTER(LEN=10) :: bas
546 : CHARACTER(LEN=default_string_length) :: name, unit_str
547 : INTEGER :: ikind, nkind, output_unit
548 : REAL(KIND=dp) :: conv, kind_radius
549 : TYPE(cp_logger_type), POINTER :: logger
550 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
551 :
552 996 : NULLIFY (logger)
553 996 : logger => cp_get_default_logger()
554 996 : NULLIFY (orb_basis_set)
555 996 : bas = "ORBITAL "
556 :
557 996 : nkind = SIZE(atomic_kind_set)
558 : ! *** Print the kind radii ***
559 : output_unit = cp_print_key_unit_nr(logger, subsys_section, &
560 996 : "PRINT%RADII/KIND_RADII", extension=".Log")
561 996 : CALL section_vals_val_get(subsys_section, "PRINT%RADII%UNIT", c_val=unit_str)
562 996 : conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
563 996 : IF (output_unit > 0) THEN
564 : WRITE (UNIT=output_unit, FMT="(/,T2,A,T56,A,T63,A,T75,A)") &
565 2 : "RADII: "//TRIM(bas)//" BASIS in "//TRIM(unit_str), &
566 4 : "Kind", "Label", "Radius"
567 8 : DO ikind = 1, nkind
568 6 : CALL get_atomic_kind(atomic_kind_set(ikind), name=name)
569 6 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
570 8 : IF (ASSOCIATED(orb_basis_set)) THEN
571 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
572 6 : kind_radius=kind_radius)
573 : WRITE (UNIT=output_unit, FMT="(T55,I5,T63,A5,T69,F12.6)") &
574 6 : ikind, name, kind_radius*conv
575 : ELSE
576 : WRITE (UNIT=output_unit, FMT="(T55,I5,T63,A5,T72,A9)") &
577 0 : ikind, name, "no basis"
578 : END IF
579 : END DO
580 : END IF
581 : CALL cp_print_key_finished_output(output_unit, logger, subsys_section, &
582 996 : "PRINT%RADII/KIND_RADII")
583 :
584 996 : END SUBROUTINE write_kind_radii
585 :
586 : ! **************************************************************************************************
587 : !> \brief Write the radii of the core charge distributions to the output
588 : !> unit.
589 : !> \param atomic_kind_set ...
590 : !> \param qs_kind_set ...
591 : !> \param subsys_section ...
592 : !> \date 15.09.2000
593 : !> \author MK
594 : !> \version 1.0
595 : ! **************************************************************************************************
596 6544 : SUBROUTINE write_core_charge_radii(atomic_kind_set, qs_kind_set, subsys_section)
597 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
598 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
599 : TYPE(section_vals_type), POINTER :: subsys_section
600 :
601 : CHARACTER(LEN=default_string_length) :: name, unit_str
602 : INTEGER :: ikind, nkind, output_unit
603 : REAL(KIND=dp) :: conv, core_charge_radius
604 : TYPE(all_potential_type), POINTER :: all_potential
605 : TYPE(cp_logger_type), POINTER :: logger
606 : TYPE(gth_potential_type), POINTER :: gth_potential
607 :
608 6544 : NULLIFY (logger)
609 6544 : logger => cp_get_default_logger()
610 : output_unit = cp_print_key_unit_nr(logger, subsys_section, &
611 6544 : "PRINT%RADII/CORE_CHARGE_RADII", extension=".Log")
612 6544 : CALL section_vals_val_get(subsys_section, "PRINT%RADII%UNIT", c_val=unit_str)
613 6544 : conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
614 6544 : IF (output_unit > 0) THEN
615 33 : nkind = SIZE(atomic_kind_set)
616 : WRITE (UNIT=output_unit, FMT="(/,T2,A,T56,A,T63,A,T75,A)") &
617 : "RADII: CORE CHARGE DISTRIBUTIONS in "// &
618 33 : TRIM(unit_str), "Kind", "Label", "Radius"
619 85 : DO ikind = 1, nkind
620 52 : CALL get_atomic_kind(atomic_kind_set(ikind), name=name)
621 : CALL get_qs_kind(qs_kind_set(ikind), &
622 52 : all_potential=all_potential, gth_potential=gth_potential)
623 :
624 85 : IF (ASSOCIATED(all_potential)) THEN
625 : CALL get_potential(potential=all_potential, &
626 22 : core_charge_radius=core_charge_radius)
627 : WRITE (UNIT=output_unit, FMT="(T55,I5,T63,A5,T69,F12.6)") &
628 22 : ikind, name, core_charge_radius*conv
629 30 : ELSE IF (ASSOCIATED(gth_potential)) THEN
630 : CALL get_potential(potential=gth_potential, &
631 30 : core_charge_radius=core_charge_radius)
632 : WRITE (UNIT=output_unit, FMT="(T55,I5,T63,A5,T69,F12.6)") &
633 30 : ikind, name, core_charge_radius*conv
634 : END IF
635 : END DO
636 : END IF
637 : CALL cp_print_key_finished_output(output_unit, logger, subsys_section, &
638 6544 : "PRINT%RADII/CORE_CHARGE_RADII")
639 6544 : END SUBROUTINE write_core_charge_radii
640 :
641 : ! **************************************************************************************************
642 : !> \brief Write the orbital basis function radii to the output unit.
643 : !> \param basis ...
644 : !> \param atomic_kind_set ...
645 : !> \param qs_kind_set ...
646 : !> \param subsys_section ...
647 : !> \date 05.06.2000
648 : !> \author MK
649 : !> \version 1.0
650 : ! **************************************************************************************************
651 19632 : SUBROUTINE write_pgf_orb_radii(basis, atomic_kind_set, qs_kind_set, subsys_section)
652 :
653 : CHARACTER(len=*) :: basis
654 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
655 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
656 : TYPE(section_vals_type), POINTER :: subsys_section
657 :
658 : CHARACTER(LEN=10) :: bas
659 : CHARACTER(LEN=default_string_length) :: name, unit_str
660 : INTEGER :: ikind, ipgf, iset, nkind, nset, &
661 : output_unit
662 19632 : INTEGER, DIMENSION(:), POINTER :: npgf
663 : REAL(KIND=dp) :: conv, kind_radius, short_kind_radius
664 19632 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius
665 19632 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pgf_radius
666 : TYPE(cp_logger_type), POINTER :: logger
667 : TYPE(gto_basis_set_type), POINTER :: aux_basis_set, lri_basis_set, &
668 : orb_basis_set
669 :
670 19632 : NULLIFY (logger)
671 39264 : logger => cp_get_default_logger()
672 19632 : NULLIFY (aux_basis_set, orb_basis_set, lri_basis_set)
673 19632 : bas = " "
674 19632 : bas(1:3) = basis(1:3)
675 19632 : CALL uppercase(bas)
676 19632 : IF (bas(1:3) == "ORB") THEN
677 6544 : bas = "ORBITAL "
678 13088 : ELSE IF (bas(1:3) == "AUX") THEN
679 6544 : bas = "AUXILLIARY"
680 6544 : ELSE IF (bas(1:3) == "LRI") THEN
681 6544 : bas = "LOCAL RI"
682 : ELSE
683 0 : CPABORT("Undefined basis set type")
684 : END IF
685 :
686 19632 : nkind = SIZE(qs_kind_set)
687 :
688 : ! *** Print the kind radii ***
689 : output_unit = cp_print_key_unit_nr(logger, subsys_section, &
690 19632 : "PRINT%RADII/KIND_RADII", extension=".Log")
691 19632 : CALL section_vals_val_get(subsys_section, "PRINT%RADII%UNIT", c_val=unit_str)
692 19632 : conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
693 19632 : IF (output_unit > 0) THEN
694 : WRITE (UNIT=output_unit, FMT="(/,T2,A,T46,A,T53,A,T63,A,T71,A)") &
695 99 : "RADII: "//TRIM(bas)//" BASIS in "//TRIM(unit_str), &
696 198 : "Kind", "Label", "Radius", "OCE Radius"
697 255 : DO ikind = 1, nkind
698 156 : CALL get_atomic_kind(atomic_kind_set(ikind), name=name)
699 156 : IF (bas(1:3) == "ORB") THEN
700 52 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
701 52 : IF (ASSOCIATED(orb_basis_set)) THEN
702 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
703 : kind_radius=kind_radius, &
704 36 : short_kind_radius=short_kind_radius)
705 : END IF
706 104 : ELSE IF (bas(1:3) == "AUX") THEN
707 52 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_basis_set, basis_type="AUX")
708 52 : IF (ASSOCIATED(aux_basis_set)) THEN
709 : CALL get_gto_basis_set(gto_basis_set=aux_basis_set, &
710 0 : kind_radius=kind_radius)
711 : END IF
712 52 : ELSE IF (bas(1:3) == "LOC") THEN
713 52 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=lri_basis_set, basis_type="LRI_AUX")
714 52 : IF (ASSOCIATED(lri_basis_set)) THEN
715 : CALL get_gto_basis_set(gto_basis_set=lri_basis_set, &
716 2 : kind_radius=kind_radius)
717 : END IF
718 : ELSE
719 0 : CPABORT("Undefined basis set type")
720 : END IF
721 255 : IF (ASSOCIATED(aux_basis_set) .OR. ASSOCIATED(orb_basis_set)) THEN
722 : WRITE (UNIT=output_unit, FMT="(T45,I5,T53,A5,T57,F12.6,T69,F12.6)") &
723 36 : ikind, name, kind_radius*conv, short_kind_radius*conv
724 : ELSE
725 : WRITE (UNIT=output_unit, FMT="(T55,I5,T63,A5,T72,A9)") &
726 120 : ikind, name, "no basis"
727 : END IF
728 : END DO
729 : END IF
730 : CALL cp_print_key_finished_output(output_unit, logger, subsys_section, &
731 19632 : "PRINT%RADII/KIND_RADII")
732 :
733 : ! *** Print the shell set radii ***
734 : output_unit = cp_print_key_unit_nr(logger, subsys_section, &
735 19632 : "PRINT%RADII/SET_RADII", extension=".Log")
736 19632 : IF (output_unit > 0) THEN
737 : WRITE (UNIT=output_unit, FMT="(/,T2,A,T51,A,T57,A,T65,A,T75,A)") &
738 : "RADII: SHELL SETS OF "//TRIM(bas)//" BASIS in "// &
739 99 : TRIM(unit_str), "Kind", "Label", "Set", "Radius"
740 255 : DO ikind = 1, nkind
741 156 : CALL get_atomic_kind(atomic_kind_set(ikind), name=name)
742 156 : IF (bas(1:3) == "ORB") THEN
743 52 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
744 52 : IF (ASSOCIATED(orb_basis_set)) THEN
745 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
746 : nset=nset, &
747 36 : set_radius=set_radius)
748 : END IF
749 104 : ELSE IF (bas(1:3) == "AUX") THEN
750 52 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_basis_set, basis_type="AUX")
751 52 : IF (ASSOCIATED(aux_basis_set)) THEN
752 : CALL get_gto_basis_set(gto_basis_set=aux_basis_set, &
753 : nset=nset, &
754 0 : set_radius=set_radius)
755 : END IF
756 52 : ELSE IF (bas(1:3) == "LOC") THEN
757 52 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=lri_basis_set, basis_type="LRI_AUX")
758 52 : IF (ASSOCIATED(lri_basis_set)) THEN
759 : CALL get_gto_basis_set(gto_basis_set=lri_basis_set, &
760 : nset=nset, &
761 2 : set_radius=set_radius)
762 : END IF
763 : ELSE
764 0 : CPABORT("Undefined basis set type")
765 : END IF
766 255 : IF (ASSOCIATED(aux_basis_set) .OR. ASSOCIATED(orb_basis_set)) THEN
767 : WRITE (UNIT=output_unit, FMT="(T50,I5,T57,A5,(T63,I5,T69,F12.6))") &
768 150 : ikind, name, (iset, set_radius(iset)*conv, iset=1, nset)
769 : ELSE
770 : WRITE (UNIT=output_unit, FMT="(T50,I5,T58,A5,T73,A8)") &
771 120 : ikind, name, "no basis"
772 : END IF
773 : END DO
774 : END IF
775 : CALL cp_print_key_finished_output(output_unit, logger, subsys_section, &
776 19632 : "PRINT%RADII/SET_RADII")
777 : ! *** Print the primitive Gaussian function radii ***
778 : output_unit = cp_print_key_unit_nr(logger, subsys_section, &
779 19632 : "PRINT%RADII/PGF_RADII", extension=".Log")
780 19632 : IF (output_unit > 0) THEN
781 : WRITE (UNIT=output_unit, FMT="(/,T2,A,T51,A,T57,A,T65,A,T75,A)") &
782 : "RADII: PRIMITIVE GAUSSIANS OF "//TRIM(bas)//" BASIS in "// &
783 99 : TRIM(unit_str), "Kind", "Label", "Set", "Radius"
784 255 : DO ikind = 1, nkind
785 156 : CALL get_atomic_kind(atomic_kind_set(ikind), name=name)
786 156 : IF (bas(1:3) == "ORB") THEN
787 52 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
788 52 : IF (ASSOCIATED(orb_basis_set)) THEN
789 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
790 : nset=nset, &
791 : npgf=npgf, &
792 36 : pgf_radius=pgf_radius)
793 : END IF
794 104 : ELSE IF (bas(1:3) == "AUX") THEN
795 52 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_basis_set, basis_type="AUX")
796 52 : IF (ASSOCIATED(aux_basis_set)) THEN
797 : CALL get_gto_basis_set(gto_basis_set=aux_basis_set, &
798 : nset=nset, &
799 : npgf=npgf, &
800 0 : pgf_radius=pgf_radius)
801 : END IF
802 52 : ELSE IF (bas(1:3) == "LOC") THEN
803 52 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=lri_basis_set, basis_type="LRI_AUX")
804 52 : IF (ASSOCIATED(lri_basis_set)) THEN
805 : CALL get_gto_basis_set(gto_basis_set=lri_basis_set, &
806 : nset=nset, &
807 : npgf=npgf, &
808 2 : pgf_radius=pgf_radius)
809 : END IF
810 : ELSE
811 0 : CPABORT("Undefined basis type")
812 : END IF
813 156 : IF (ASSOCIATED(aux_basis_set) .OR. ASSOCIATED(orb_basis_set) .OR. &
814 99 : ASSOCIATED(lri_basis_set)) THEN
815 177 : DO iset = 1, nset
816 : WRITE (UNIT=output_unit, FMT="(T50,I5,T57,A5,T63,I5,(T69,F12.6))") &
817 139 : ikind, name, iset, &
818 586 : (pgf_radius(ipgf, iset)*conv, ipgf=1, npgf(iset))
819 : END DO
820 : ELSE
821 : WRITE (UNIT=output_unit, FMT="(T50,I5,T58,A5,T73,A8)") &
822 118 : ikind, name, "no basis"
823 : END IF
824 : END DO
825 : END IF
826 : CALL cp_print_key_finished_output(output_unit, logger, subsys_section, &
827 19632 : "PRINT%RADII/PGF_RADII")
828 19632 : END SUBROUTINE write_pgf_orb_radii
829 :
830 : ! **************************************************************************************************
831 : !> \brief Write the radii of the exponential functions of the Goedecker
832 : !> pseudopotential (GTH, local part) to the logical unit number
833 : !> "output_unit".
834 : !> \param atomic_kind_set ...
835 : !> \param qs_kind_set ...
836 : !> \param subsys_section ...
837 : !> \date 06.10.2000
838 : !> \author MK
839 : !> \version 1.0
840 : ! **************************************************************************************************
841 6544 : SUBROUTINE write_ppl_radii(atomic_kind_set, qs_kind_set, subsys_section)
842 :
843 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
844 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
845 : TYPE(section_vals_type), POINTER :: subsys_section
846 :
847 : CHARACTER(LEN=default_string_length) :: name, unit_str
848 : INTEGER :: ikind, nkind, output_unit
849 : REAL(KIND=dp) :: conv, ppl_radius
850 : TYPE(cp_logger_type), POINTER :: logger
851 : TYPE(gth_potential_type), POINTER :: gth_potential
852 :
853 6544 : NULLIFY (logger)
854 6544 : logger => cp_get_default_logger()
855 : output_unit = cp_print_key_unit_nr(logger, subsys_section, &
856 6544 : "PRINT%RADII/PPL_RADII", extension=".Log")
857 6544 : CALL section_vals_val_get(subsys_section, "PRINT%RADII%UNIT", c_val=unit_str)
858 6544 : conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
859 6544 : IF (output_unit > 0) THEN
860 33 : nkind = SIZE(atomic_kind_set)
861 : WRITE (UNIT=output_unit, FMT="(/,T2,A,T56,A,T63,A,T75,A)") &
862 : "RADII: LOCAL PART OF GTH/ELP PP in "// &
863 33 : TRIM(unit_str), "Kind", "Label", "Radius"
864 85 : DO ikind = 1, nkind
865 52 : CALL get_atomic_kind(atomic_kind_set(ikind), name=name)
866 52 : CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential)
867 85 : IF (ASSOCIATED(gth_potential)) THEN
868 : CALL get_potential(potential=gth_potential, &
869 30 : ppl_radius=ppl_radius)
870 : WRITE (UNIT=output_unit, FMT="(T55,I5,T63,A5,T69,F12.6)") &
871 30 : ikind, name, ppl_radius*conv
872 : END IF
873 : END DO
874 : END IF
875 : CALL cp_print_key_finished_output(output_unit, logger, subsys_section, &
876 6544 : "PRINT%RADII/PPL_RADII")
877 6544 : END SUBROUTINE write_ppl_radii
878 :
879 : ! **************************************************************************************************
880 : !> \brief Write the radii of the projector functions of the Goedecker
881 : !> pseudopotential (GTH, non-local part) to the logical unit number
882 : !> "output_unit".
883 : !> \param atomic_kind_set ...
884 : !> \param qs_kind_set ...
885 : !> \param subsys_section ...
886 : !> \date 06.10.2000
887 : !> \author MK
888 : !> \version 1.0
889 : ! **************************************************************************************************
890 6544 : SUBROUTINE write_ppnl_radii(atomic_kind_set, qs_kind_set, subsys_section)
891 :
892 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
893 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
894 : TYPE(section_vals_type), POINTER :: subsys_section
895 :
896 : CHARACTER(LEN=default_string_length) :: name, unit_str
897 : INTEGER :: ikind, nkind, output_unit
898 : REAL(KIND=dp) :: conv, ppnl_radius
899 : TYPE(cp_logger_type), POINTER :: logger
900 : TYPE(gth_potential_type), POINTER :: gth_potential
901 :
902 6544 : NULLIFY (logger)
903 6544 : logger => cp_get_default_logger()
904 : output_unit = cp_print_key_unit_nr(logger, subsys_section, &
905 6544 : "PRINT%RADII/PPNL_RADII", extension=".Log")
906 6544 : CALL section_vals_val_get(subsys_section, "PRINT%RADII%UNIT", c_val=unit_str)
907 6544 : conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
908 6544 : IF (output_unit > 0) THEN
909 33 : nkind = SIZE(atomic_kind_set)
910 : WRITE (UNIT=output_unit, FMT="(/,T2,A,T56,A,T63,A,T75,A)") &
911 : "RADII: NON-LOCAL PART OF GTH PP in "// &
912 33 : TRIM(unit_str), "Kind", "Label", "Radius"
913 85 : DO ikind = 1, nkind
914 52 : CALL get_atomic_kind(atomic_kind_set(ikind), name=name)
915 52 : CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential)
916 85 : IF (ASSOCIATED(gth_potential)) THEN
917 : CALL get_potential(potential=gth_potential, &
918 30 : ppnl_radius=ppnl_radius)
919 : WRITE (UNIT=output_unit, FMT="(T55,I5,T63,A5,T69,F12.6)") &
920 30 : ikind, name, ppnl_radius*conv
921 : END IF
922 : END DO
923 : END IF
924 : CALL cp_print_key_finished_output(output_unit, logger, subsys_section, &
925 6544 : "PRINT%RADII/PPNL_RADII")
926 6544 : END SUBROUTINE write_ppnl_radii
927 :
928 : ! **************************************************************************************************
929 : !> \brief Write the radii of the one center projector
930 : !> \param atomic_kind_set ...
931 : !> \param qs_kind_set ...
932 : !> \param subsys_section ...
933 : !> \version 1.0
934 : ! **************************************************************************************************
935 6544 : SUBROUTINE write_paw_radii(atomic_kind_set, qs_kind_set, subsys_section)
936 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
937 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
938 : TYPE(section_vals_type), POINTER :: subsys_section
939 :
940 : CHARACTER(LEN=default_string_length) :: name, unit_str
941 : INTEGER :: ikind, nkind, output_unit
942 : LOGICAL :: paw_atom
943 : REAL(KIND=dp) :: conv, rcprj
944 : TYPE(cp_logger_type), POINTER :: logger
945 : TYPE(paw_proj_set_type), POINTER :: paw_proj_set
946 :
947 6544 : NULLIFY (logger)
948 6544 : logger => cp_get_default_logger()
949 : output_unit = cp_print_key_unit_nr(logger, subsys_section, &
950 6544 : "PRINT%RADII/GAPW_PRJ_RADII", extension=".Log")
951 6544 : CALL section_vals_val_get(subsys_section, "PRINT%RADII%UNIT", c_val=unit_str)
952 6544 : conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
953 6544 : IF (output_unit > 0) THEN
954 33 : nkind = SIZE(qs_kind_set)
955 : WRITE (UNIT=output_unit, FMT="(/,T2,A,T56,A,T63,A,T75,A)") &
956 : "RADII: ONE CENTER PROJECTORS in "// &
957 33 : TRIM(unit_str), "Kind", "Label", "Radius"
958 85 : DO ikind = 1, nkind
959 52 : CALL get_atomic_kind(atomic_kind_set(ikind), name=name)
960 : CALL get_qs_kind(qs_kind_set(ikind), &
961 52 : paw_atom=paw_atom, paw_proj_set=paw_proj_set)
962 85 : IF (paw_atom .AND. ASSOCIATED(paw_proj_set)) THEN
963 : CALL get_paw_proj_set(paw_proj_set=paw_proj_set, &
964 0 : rcprj=rcprj)
965 : WRITE (UNIT=output_unit, FMT="(T55,I5,T63,A5,T69,F12.6)") &
966 0 : ikind, name, rcprj*conv
967 : END IF
968 : END DO
969 : END IF
970 : CALL cp_print_key_finished_output(output_unit, logger, subsys_section, &
971 6544 : "PRINT%RADII/GAPW_PRJ_RADII")
972 6544 : END SUBROUTINE write_paw_radii
973 :
974 : END MODULE qs_interactions
|
