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