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