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