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 Define the quickstep kind type and their sub types
10 : !> \author Ole Schuett
11 : !>
12 : !> <b>Modification history:</b>
13 : !> - 01.2002 creation [MK]
14 : !> - 04.2002 added pao [fawzi]
15 : !> - 09.2002 adapted for POL/KG use [GT]
16 : !> - 02.2004 flexible normalization of basis sets [jgh]
17 : !> - 03.2004 attach/detach routines [jgh]
18 : !> - 10.2004 removed pao [fawzi]
19 : !> - 08.2014 separated qs-related stuff from atomic_kind_types.F [Ole Schuett]
20 : !> - 07.2015 new container for basis sets [jgh]
21 : !> - 04.2021 init dft_plus_u_type [MK]
22 : ! **************************************************************************************************
23 : MODULE qs_kind_types
24 : USE atom_sgp, ONLY: atom_sgp_potential_type,&
25 : atom_sgp_release,&
26 : sgp_construction
27 : USE atom_types, ONLY: atom_ecppot_type,&
28 : lmat,&
29 : read_ecp_potential
30 : USE atom_upf, ONLY: atom_read_upf,&
31 : atom_release_upf,&
32 : atom_upfpot_type
33 : USE atomic_kind_types, ONLY: atomic_kind_type,&
34 : get_atomic_kind
35 : USE basis_set_container_types, ONLY: add_basis_set_to_container,&
36 : basis_set_container_type,&
37 : get_basis_from_container,&
38 : remove_basis_from_container,&
39 : remove_basis_set_container
40 : USE basis_set_types, ONLY: &
41 : allocate_gto_basis_set, allocate_sto_basis_set, combine_basis_sets, &
42 : create_gto_from_sto_basis, deallocate_sto_basis_set, get_gto_basis_set, &
43 : gto_basis_set_type, init_aux_basis_set, init_orb_basis_set, read_gto_basis_set, &
44 : read_sto_basis_set, sto_basis_set_type, write_gto_basis_set, write_orb_basis_set
45 : USE cp_control_types, ONLY: dft_control_type,&
46 : qs_control_type,&
47 : xtb_control_type
48 : USE cp_log_handling, ONLY: cp_get_default_logger,&
49 : cp_logger_get_default_io_unit,&
50 : cp_logger_type
51 : USE cp_output_handling, ONLY: cp_p_file,&
52 : cp_print_key_finished_output,&
53 : cp_print_key_should_output,&
54 : cp_print_key_unit_nr
55 : USE external_potential_types, ONLY: &
56 : all_potential_type, allocate_potential, deallocate_potential, get_potential, &
57 : gth_potential_type, init_potential, local_potential_type, read_potential, &
58 : set_default_all_potential, set_potential, sgp_potential_type, write_potential
59 : USE gapw_1c_basis_set, ONLY: create_1c_basis
60 : USE input_constants, ONLY: &
61 : do_method_am1, do_method_dftb, do_method_mndo, do_method_mndod, do_method_pdg, &
62 : do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_pnnl, do_method_pw, &
63 : do_method_rm1, do_method_xtb, do_qs, do_sirius, gapw_1c_large, gapw_1c_medium, &
64 : gapw_1c_orb, gapw_1c_small, gapw_1c_very_large
65 : USE input_section_types, ONLY: section_vals_get,&
66 : section_vals_get_subs_vals,&
67 : section_vals_type,&
68 : section_vals_val_get
69 : USE kinds, ONLY: default_path_length,&
70 : default_string_length,&
71 : dp
72 : USE mathconstants, ONLY: pi
73 : USE message_passing, ONLY: mp_para_env_type
74 : USE orbital_pointers, ONLY: init_orbital_pointers,&
75 : nco,&
76 : ncoset
77 : USE paw_proj_set_types, ONLY: allocate_paw_proj_set,&
78 : deallocate_paw_proj_set,&
79 : get_paw_proj_set,&
80 : paw_proj_set_type,&
81 : projectors
82 : USE periodic_table, ONLY: get_ptable_info,&
83 : ptable
84 : USE physcon, ONLY: angstrom,&
85 : bohr,&
86 : evolt
87 : USE qs_dftb_types, ONLY: qs_dftb_atom_type
88 : USE qs_dftb_utils, ONLY: deallocate_dftb_atom_param,&
89 : get_dftb_atom_param,&
90 : write_dftb_atom_param
91 : USE qs_dispersion_types, ONLY: qs_atom_dispersion_type
92 : USE qs_grid_atom, ONLY: allocate_grid_atom,&
93 : deallocate_grid_atom,&
94 : grid_atom_type
95 : USE qs_harmonics_atom, ONLY: allocate_harmonics_atom,&
96 : deallocate_harmonics_atom,&
97 : harmonics_atom_type
98 : USE semi_empirical_types, ONLY: get_se_param,&
99 : semi_empirical_create,&
100 : semi_empirical_release,&
101 : semi_empirical_type,&
102 : write_se_param
103 : USE semi_empirical_utils, ONLY: init_se_param,&
104 : se_param_set_default
105 : USE soft_basis_set, ONLY: create_soft_basis
106 : USE string_utilities, ONLY: uppercase
107 : USE xtb_parameters, ONLY: xtb_set_kab
108 : USE xtb_types, ONLY: deallocate_xtb_atom_param,&
109 : get_xtb_atom_param,&
110 : write_xtb_atom_param,&
111 : xtb_atom_type
112 : #include "./base/base_uses.f90"
113 :
114 : IMPLICIT NONE
115 :
116 : PRIVATE
117 :
118 : ! Global parameters (only in this module)
119 :
120 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_kind_types'
121 :
122 : ! **************************************************************************************************
123 : !> \brief Input parameters for the DFT+U method
124 : ! **************************************************************************************************
125 : TYPE dft_plus_u_type
126 : INTEGER :: l = -1
127 : INTEGER :: n = -1
128 : INTEGER :: max_scf = -1
129 : REAL(KIND=dp) :: eps_u_ramping = 0.0_dp
130 : REAL(KIND=dp) :: eps_scf = HUGE(0.0_dp)
131 : REAL(KIND=dp) :: u_minus_j_target = 0.0_dp
132 : REAL(KIND=dp) :: u_minus_j = 0.0_dp
133 : REAL(KIND=dp) :: u_ramping = 0.0_dp
134 : REAL(KIND=dp) :: U = 0.0_dp
135 : REAL(KIND=dp) :: J = 0.0_dp
136 : REAL(KIND=dp) :: alpha = 0.0_dp
137 : REAL(KIND=dp) :: beta = 0.0_dp
138 : REAL(KIND=dp) :: J0 = 0.0_dp
139 : REAL(KIND=dp) :: occupation = -1.0_dp
140 : INTEGER, DIMENSION(:), POINTER :: orbitals => Null()
141 : LOGICAL :: init_u_ramping_each_scf = .FALSE.
142 : LOGICAL :: smear = .FALSE.
143 : REAL(KIND=dp), DIMENSION(:), POINTER :: nelec => Null()
144 : END TYPE dft_plus_u_type
145 :
146 : ! **************************************************************************************************
147 : !> \brief Holds information about a PAO potential
148 : ! **************************************************************************************************
149 : TYPE pao_potential_type
150 : INTEGER :: maxl = -1
151 : REAL(KIND=dp) :: beta = 0.0_dp
152 : REAL(KIND=dp) :: weight = 0.0_dp
153 : INTEGER :: max_projector = -1
154 : REAL(KIND=dp) :: beta_radius = HUGE(dp)
155 : END TYPE pao_potential_type
156 :
157 : ! **************************************************************************************************
158 : !> \brief Holds information about a PAO descriptor
159 : ! **************************************************************************************************
160 : TYPE pao_descriptor_type
161 : REAL(KIND=dp) :: beta = 0.0_dp
162 : REAL(KIND=dp) :: beta_radius = HUGE(dp)
163 : REAL(KIND=dp) :: weight = 0.0_dp
164 : REAL(KIND=dp) :: screening = 0.0_dp
165 : REAL(KIND=dp) :: screening_radius = HUGE(dp)
166 : END TYPE pao_descriptor_type
167 :
168 : ! **************************************************************************************************
169 : !> \brief Provides all information about a quickstep kind
170 : ! **************************************************************************************************
171 : TYPE qs_kind_type
172 : CHARACTER(LEN=default_string_length) :: name = ""
173 : CHARACTER(LEN=2) :: element_symbol = ""
174 : INTEGER :: natom = -1
175 : TYPE(all_potential_type), POINTER :: all_potential => Null()
176 : TYPE(local_potential_type), POINTER :: tnadd_potential => Null()
177 : TYPE(gth_potential_type), POINTER :: gth_potential => Null()
178 : TYPE(sgp_potential_type), POINTER :: sgp_potential => Null()
179 : TYPE(semi_empirical_type), POINTER :: se_parameter => Null()
180 : TYPE(qs_dftb_atom_type), POINTER :: dftb_parameter => Null()
181 : TYPE(xtb_atom_type), POINTER :: xtb_parameter => Null()
182 : !
183 : TYPE(atom_upfpot_type), POINTER :: upf_potential => Null()
184 : !
185 : TYPE(basis_set_container_type), &
186 : DIMENSION(20) :: basis_sets = basis_set_container_type()
187 : ! Atomic radii
188 : REAL(KIND=dp) :: covalent_radius = 0.0_dp
189 : REAL(KIND=dp) :: vdw_radius = 0.0_dp
190 : ! GAPW specific data
191 : TYPE(paw_proj_set_type), POINTER :: paw_proj_set => Null()
192 : REAL(KIND=dp) :: hard_radius = 0.8_dp*bohr ! for hard and soft exp
193 : REAL(KIND=dp) :: hard0_radius = 0.8_dp*bohr ! for hard exp of rho0
194 : REAL(KIND=dp) :: max_rad_local = 13.2_dp*bohr ! max GTO radius used in GAPW
195 : LOGICAL :: paw_atom = .FALSE. ! needs atomic rho1
196 : LOGICAL :: gpw_type_forced = .FALSE. ! gpw atom even if with hard exponents
197 : !
198 : LOGICAL :: ghost = .FALSE.
199 : LOGICAL :: floating = .FALSE.
200 : INTEGER :: lmax_dftb = -1
201 : REAL(KIND=dp) :: dudq_dftb3 = 0.0_dp
202 : REAL(KIND=dp) :: magnetization = 0.0_dp
203 : INTEGER, DIMENSION(:, :), POINTER :: addel => Null()
204 : INTEGER, DIMENSION(:, :), POINTER :: laddel => Null()
205 : INTEGER, DIMENSION(:, :), POINTER :: naddel => Null()
206 : TYPE(harmonics_atom_type), POINTER :: harmonics => Null()
207 : TYPE(grid_atom_type), POINTER :: grid_atom => Null()
208 : INTEGER :: ngrid_rad = 50
209 : INTEGER :: ngrid_ang = 50
210 : INTEGER :: lmax_rho0 = 0
211 : INTEGER :: mao = -1
212 : INTEGER, DIMENSION(:), POINTER :: elec_conf => Null() ! used to set up the initial atomic guess
213 : LOGICAL :: bs_occupation = .FALSE.
214 : TYPE(dft_plus_u_type), POINTER :: dft_plus_u => Null()
215 : LOGICAL :: no_optimize = .TRUE.
216 : !
217 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: nlcc_pot => Null()
218 : !
219 : TYPE(qs_atom_dispersion_type), POINTER :: dispersion => Null()
220 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: reltmat => Null()
221 : INTEGER :: pao_basis_size = -1
222 : CHARACTER(LEN=default_path_length) :: pao_model_file = ""
223 : TYPE(pao_potential_type), DIMENSION(:), POINTER :: pao_potentials => Null()
224 : TYPE(pao_descriptor_type), DIMENSION(:), POINTER :: pao_descriptors => Null()
225 : END TYPE qs_kind_type
226 :
227 : ! **************************************************************************************************
228 : !> \brief Provides a vector of pointers of type qs_kind_type
229 : ! **************************************************************************************************
230 : TYPE qs_kind_p_type
231 : TYPE(qs_kind_type), DIMENSION(:), &
232 : POINTER :: qs_kind_set => NULL()
233 : END TYPE qs_kind_p_type
234 :
235 : ! Public subroutines
236 :
237 : PUBLIC :: check_qs_kind_set, &
238 : deallocate_qs_kind_set, &
239 : get_qs_kind, &
240 : get_qs_kind_set, &
241 : has_nlcc, &
242 : init_qs_kind_set, &
243 : init_gapw_basis_set, &
244 : init_gapw_nlcc, &
245 : create_qs_kind_set, &
246 : set_qs_kind, &
247 : write_qs_kind_set, &
248 : write_gto_basis_sets, &
249 : init_atom_electronic_state, set_pseudo_state
250 :
251 : ! Public data types
252 : PUBLIC :: qs_kind_type, pao_potential_type, pao_descriptor_type
253 :
254 : CONTAINS
255 :
256 : ! **************************************************************************************************
257 : !> \brief Destructor routine for a set of qs kinds
258 : !> \param qs_kind_set ...
259 : !> \date 02.01.2002
260 : !> \author Matthias Krack (MK)
261 : !> \version 2.0
262 : ! **************************************************************************************************
263 6814 : SUBROUTINE deallocate_qs_kind_set(qs_kind_set)
264 :
265 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
266 :
267 : INTEGER :: ikind, nkind
268 :
269 6814 : IF (ASSOCIATED(qs_kind_set)) THEN
270 :
271 6814 : nkind = SIZE(qs_kind_set)
272 :
273 19799 : DO ikind = 1, nkind
274 12985 : IF (ASSOCIATED(qs_kind_set(ikind)%all_potential)) THEN
275 4466 : CALL deallocate_potential(qs_kind_set(ikind)%all_potential)
276 : END IF
277 12985 : IF (ASSOCIATED(qs_kind_set(ikind)%tnadd_potential)) THEN
278 20 : CALL deallocate_potential(qs_kind_set(ikind)%tnadd_potential)
279 : END IF
280 12985 : IF (ASSOCIATED(qs_kind_set(ikind)%gth_potential)) THEN
281 8331 : CALL deallocate_potential(qs_kind_set(ikind)%gth_potential)
282 : END IF
283 12985 : IF (ASSOCIATED(qs_kind_set(ikind)%sgp_potential)) THEN
284 24 : CALL deallocate_potential(qs_kind_set(ikind)%sgp_potential)
285 : END IF
286 12985 : IF (ASSOCIATED(qs_kind_set(ikind)%upf_potential)) THEN
287 20 : CALL atom_release_upf(qs_kind_set(ikind)%upf_potential)
288 20 : DEALLOCATE (qs_kind_set(ikind)%upf_potential)
289 : END IF
290 12985 : IF (ASSOCIATED(qs_kind_set(ikind)%se_parameter)) THEN
291 2240 : CALL semi_empirical_release(qs_kind_set(ikind)%se_parameter)
292 : END IF
293 12985 : IF (ASSOCIATED(qs_kind_set(ikind)%dftb_parameter)) THEN
294 480 : CALL deallocate_dftb_atom_param(qs_kind_set(ikind)%dftb_parameter)
295 : END IF
296 12985 : IF (ASSOCIATED(qs_kind_set(ikind)%xtb_parameter)) THEN
297 692 : CALL deallocate_xtb_atom_param(qs_kind_set(ikind)%xtb_parameter)
298 : END IF
299 12985 : IF (ASSOCIATED(qs_kind_set(ikind)%paw_proj_set)) THEN
300 1618 : CALL deallocate_paw_proj_set(qs_kind_set(ikind)%paw_proj_set)
301 : END IF
302 12985 : IF (ASSOCIATED(qs_kind_set(ikind)%harmonics)) THEN
303 1928 : CALL deallocate_harmonics_atom(qs_kind_set(ikind)%harmonics)
304 : END IF
305 12985 : IF (ASSOCIATED(qs_kind_set(ikind)%grid_atom)) THEN
306 1928 : CALL deallocate_grid_atom(qs_kind_set(ikind)%grid_atom)
307 : END IF
308 12985 : IF (ASSOCIATED(qs_kind_set(ikind)%elec_conf)) THEN
309 12681 : DEALLOCATE (qs_kind_set(ikind)%elec_conf)
310 : END IF
311 :
312 12985 : IF (ASSOCIATED(qs_kind_set(ikind)%dft_plus_u)) THEN
313 32 : IF (ASSOCIATED(qs_kind_set(ikind)%dft_plus_u%orbitals)) THEN
314 4 : DEALLOCATE (qs_kind_set(ikind)%dft_plus_u%orbitals)
315 : END IF
316 32 : IF (ASSOCIATED(qs_kind_set(ikind)%dft_plus_u%nelec)) THEN
317 4 : DEALLOCATE (qs_kind_set(ikind)%dft_plus_u%nelec)
318 : END IF
319 32 : DEALLOCATE (qs_kind_set(ikind)%dft_plus_u)
320 : END IF
321 :
322 12985 : IF (ASSOCIATED(qs_kind_set(ikind)%nlcc_pot)) THEN
323 2 : DEALLOCATE (qs_kind_set(ikind)%nlcc_pot)
324 : END IF
325 :
326 12985 : IF (ASSOCIATED(qs_kind_set(ikind)%dispersion)) THEN
327 902 : DEALLOCATE (qs_kind_set(ikind)%dispersion)
328 : END IF
329 12985 : IF (ASSOCIATED(qs_kind_set(ikind)%addel)) THEN
330 60 : DEALLOCATE (qs_kind_set(ikind)%addel)
331 : END IF
332 12985 : IF (ASSOCIATED(qs_kind_set(ikind)%naddel)) THEN
333 60 : DEALLOCATE (qs_kind_set(ikind)%naddel)
334 : END IF
335 12985 : IF (ASSOCIATED(qs_kind_set(ikind)%laddel)) THEN
336 60 : DEALLOCATE (qs_kind_set(ikind)%laddel)
337 : END IF
338 12985 : IF (ASSOCIATED(qs_kind_set(ikind)%reltmat)) THEN
339 26 : DEALLOCATE (qs_kind_set(ikind)%reltmat)
340 : END IF
341 :
342 12985 : IF (ASSOCIATED(qs_kind_set(ikind)%pao_potentials)) THEN
343 9571 : DEALLOCATE (qs_kind_set(ikind)%pao_potentials)
344 : END IF
345 12985 : IF (ASSOCIATED(qs_kind_set(ikind)%pao_descriptors)) THEN
346 9571 : DEALLOCATE (qs_kind_set(ikind)%pao_descriptors)
347 : END IF
348 :
349 19799 : CALL remove_basis_set_container(qs_kind_set(ikind)%basis_sets)
350 :
351 : END DO
352 6814 : DEALLOCATE (qs_kind_set)
353 : ELSE
354 : CALL cp_abort(__LOCATION__, &
355 : "The pointer qs_kind_set is not associated and "// &
356 0 : "cannot be deallocated")
357 : END IF
358 :
359 6814 : END SUBROUTINE deallocate_qs_kind_set
360 :
361 : ! **************************************************************************************************
362 : !> \brief Get attributes of an atomic kind.
363 : !> \param qs_kind ...
364 : !> \param basis_set ...
365 : !> \param basis_type ...
366 : !> \param ncgf ...
367 : !> \param nsgf ...
368 : !> \param all_potential ...
369 : !> \param tnadd_potential ...
370 : !> \param gth_potential ...
371 : !> \param sgp_potential ...
372 : !> \param upf_potential ...
373 : !> \param se_parameter ...
374 : !> \param dftb_parameter ...
375 : !> \param xtb_parameter ...
376 : !> \param dftb3_param ...
377 : !> \param zatom ...
378 : !> \param zeff ...
379 : !> \param elec_conf ...
380 : !> \param mao ...
381 : !> \param lmax_dftb ...
382 : !> \param alpha_core_charge ...
383 : !> \param ccore_charge ...
384 : !> \param core_charge ...
385 : !> \param core_charge_radius ...
386 : !> \param paw_proj_set ...
387 : !> \param paw_atom ...
388 : !> \param hard_radius ...
389 : !> \param hard0_radius ...
390 : !> \param max_rad_local ...
391 : !> \param covalent_radius ...
392 : !> \param vdw_radius ...
393 : !> \param gpw_type_forced ...
394 : !> \param harmonics ...
395 : !> \param max_iso_not0 ...
396 : !> \param max_s_harm ...
397 : !> \param grid_atom ...
398 : !> \param ngrid_ang ...
399 : !> \param ngrid_rad ...
400 : !> \param lmax_rho0 ...
401 : !> \param dft_plus_u_atom ...
402 : !> \param l_of_dft_plus_u ...
403 : !> \param n_of_dft_plus_u ...
404 : !> \param u_minus_j ...
405 : !> \param U_of_dft_plus_u ...
406 : !> \param J_of_dft_plus_u ...
407 : !> \param alpha_of_dft_plus_u ...
408 : !> \param beta_of_dft_plus_u ...
409 : !> \param J0_of_dft_plus_u ...
410 : !> \param occupation_of_dft_plus_u ...
411 : !> \param dispersion ...
412 : !> \param bs_occupation ...
413 : !> \param magnetization ...
414 : !> \param no_optimize ...
415 : !> \param addel ...
416 : !> \param laddel ...
417 : !> \param naddel ...
418 : !> \param orbitals ...
419 : !> \param max_scf ...
420 : !> \param eps_scf ...
421 : !> \param smear ...
422 : !> \param u_ramping ...
423 : !> \param u_minus_j_target ...
424 : !> \param eps_u_ramping ...
425 : !> \param init_u_ramping_each_scf ...
426 : !> \param reltmat ...
427 : !> \param ghost ...
428 : !> \param floating ...
429 : !> \param name ...
430 : !> \param element_symbol ...
431 : !> \param pao_basis_size ...
432 : !> \param pao_model_file ...
433 : !> \param pao_potentials ...
434 : !> \param pao_descriptors ...
435 : !> \param nelec ...
436 : ! **************************************************************************************************
437 53059330 : SUBROUTINE get_qs_kind(qs_kind, &
438 : basis_set, basis_type, ncgf, nsgf, &
439 : all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, &
440 : se_parameter, dftb_parameter, xtb_parameter, &
441 : dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, &
442 : alpha_core_charge, ccore_charge, core_charge, core_charge_radius, &
443 : paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, &
444 : covalent_radius, vdw_radius, &
445 : gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, &
446 : ngrid_ang, ngrid_rad, lmax_rho0, &
447 : dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, &
448 : u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, &
449 : alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, &
450 : bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, &
451 : max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, &
452 : init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, &
453 : pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
454 :
455 : TYPE(qs_kind_type) :: qs_kind
456 : TYPE(gto_basis_set_type), OPTIONAL, POINTER :: basis_set
457 : CHARACTER(len=*), OPTIONAL :: basis_type
458 : INTEGER, INTENT(OUT), OPTIONAL :: ncgf, nsgf
459 : TYPE(all_potential_type), OPTIONAL, POINTER :: all_potential
460 : TYPE(local_potential_type), OPTIONAL, POINTER :: tnadd_potential
461 : TYPE(gth_potential_type), OPTIONAL, POINTER :: gth_potential
462 : TYPE(sgp_potential_type), OPTIONAL, POINTER :: sgp_potential
463 : TYPE(atom_upfpot_type), OPTIONAL, POINTER :: upf_potential
464 : TYPE(semi_empirical_type), OPTIONAL, POINTER :: se_parameter
465 : TYPE(qs_dftb_atom_type), OPTIONAL, POINTER :: dftb_parameter
466 : TYPE(xtb_atom_type), OPTIONAL, POINTER :: xtb_parameter
467 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: dftb3_param
468 : INTEGER, INTENT(OUT), OPTIONAL :: zatom
469 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: zeff
470 : INTEGER, DIMENSION(:), OPTIONAL, POINTER :: elec_conf
471 : INTEGER, INTENT(OUT), OPTIONAL :: mao, lmax_dftb
472 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: alpha_core_charge, ccore_charge, &
473 : core_charge, core_charge_radius
474 : TYPE(paw_proj_set_type), OPTIONAL, POINTER :: paw_proj_set
475 : LOGICAL, INTENT(OUT), OPTIONAL :: paw_atom
476 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: hard_radius, hard0_radius, &
477 : max_rad_local, covalent_radius, &
478 : vdw_radius
479 : LOGICAL, INTENT(OUT), OPTIONAL :: gpw_type_forced
480 : TYPE(harmonics_atom_type), OPTIONAL, POINTER :: harmonics
481 : INTEGER, INTENT(OUT), OPTIONAL :: max_iso_not0, max_s_harm
482 : TYPE(grid_atom_type), OPTIONAL, POINTER :: grid_atom
483 : INTEGER, INTENT(OUT), OPTIONAL :: ngrid_ang, ngrid_rad, lmax_rho0
484 : LOGICAL, INTENT(OUT), OPTIONAL :: dft_plus_u_atom
485 : INTEGER, INTENT(OUT), OPTIONAL :: l_of_dft_plus_u, n_of_dft_plus_u
486 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, &
487 : alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u
488 : TYPE(qs_atom_dispersion_type), OPTIONAL, POINTER :: dispersion
489 : LOGICAL, INTENT(OUT), OPTIONAL :: bs_occupation
490 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: magnetization
491 : LOGICAL, INTENT(OUT), OPTIONAL :: no_optimize
492 : INTEGER, DIMENSION(:, :), OPTIONAL, POINTER :: addel, laddel, naddel
493 : INTEGER, DIMENSION(:), OPTIONAL, POINTER :: orbitals
494 : INTEGER, OPTIONAL :: max_scf
495 : REAL(KIND=dp), OPTIONAL :: eps_scf
496 : LOGICAL, OPTIONAL :: smear
497 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: u_ramping, u_minus_j_target, &
498 : eps_u_ramping
499 : LOGICAL, OPTIONAL :: init_u_ramping_each_scf
500 : REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: reltmat
501 : LOGICAL, OPTIONAL :: ghost, floating
502 : CHARACTER(LEN=default_string_length), &
503 : INTENT(OUT), OPTIONAL :: name
504 : CHARACTER(LEN=2), INTENT(OUT), OPTIONAL :: element_symbol
505 : INTEGER, INTENT(OUT), OPTIONAL :: pao_basis_size
506 : CHARACTER(LEN=default_path_length), INTENT(OUT), &
507 : OPTIONAL :: pao_model_file
508 : TYPE(pao_potential_type), DIMENSION(:), OPTIONAL, &
509 : POINTER :: pao_potentials
510 : TYPE(pao_descriptor_type), DIMENSION(:), &
511 : OPTIONAL, POINTER :: pao_descriptors
512 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: nelec
513 :
514 : CHARACTER(LEN=default_string_length) :: my_basis_type
515 : INTEGER :: l
516 : LOGICAL :: found
517 : TYPE(gto_basis_set_type), POINTER :: tmp_basis_set
518 :
519 : ! Retrieve basis set from the kind container
520 53059330 : IF (PRESENT(basis_type)) THEN
521 7520976 : my_basis_type = basis_type
522 : ELSE
523 45538354 : my_basis_type = "ORB"
524 : END IF
525 :
526 53059330 : IF (PRESENT(basis_set)) THEN
527 : CALL get_basis_from_container(qs_kind%basis_sets, basis_set=basis_set, &
528 8131327 : basis_type=my_basis_type)
529 : END IF
530 :
531 53059330 : IF (PRESENT(ncgf)) THEN
532 : CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis_set, &
533 960 : basis_type=my_basis_type)
534 960 : IF (ASSOCIATED(tmp_basis_set)) THEN
535 960 : CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, ncgf=ncgf)
536 0 : ELSE IF (ASSOCIATED(qs_kind%dftb_parameter)) THEN
537 0 : l = qs_kind%dftb_parameter%lmax
538 0 : ncgf = ((l + 1)*(l + 2)*(l + 3))/6
539 : ELSE
540 0 : ncgf = 0
541 : END IF
542 : END IF
543 :
544 53059330 : IF (PRESENT(nsgf)) THEN
545 : CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis_set, &
546 262629 : basis_type=my_basis_type)
547 262629 : IF (ASSOCIATED(tmp_basis_set)) THEN
548 160815 : CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, nsgf=nsgf)
549 101814 : ELSE IF (ASSOCIATED(qs_kind%dftb_parameter)) THEN
550 101810 : nsgf = qs_kind%dftb_parameter%natorb
551 : ELSE
552 4 : nsgf = 0
553 : END IF
554 : END IF
555 :
556 53059330 : IF (PRESENT(all_potential)) all_potential => qs_kind%all_potential
557 53059330 : IF (PRESENT(tnadd_potential)) tnadd_potential => qs_kind%tnadd_potential
558 53059330 : IF (PRESENT(gth_potential)) gth_potential => qs_kind%gth_potential
559 53059330 : IF (PRESENT(sgp_potential)) sgp_potential => qs_kind%sgp_potential
560 53059330 : IF (PRESENT(upf_potential)) upf_potential => qs_kind%upf_potential
561 53059330 : IF (PRESENT(se_parameter)) se_parameter => qs_kind%se_parameter
562 53059330 : IF (PRESENT(dftb_parameter)) dftb_parameter => qs_kind%dftb_parameter
563 53059330 : IF (PRESENT(xtb_parameter)) xtb_parameter => qs_kind%xtb_parameter
564 :
565 53059330 : IF (PRESENT(element_symbol)) element_symbol = qs_kind%element_symbol
566 53059330 : IF (PRESENT(name)) name = qs_kind%name
567 53059330 : IF (PRESENT(dftb3_param)) dftb3_param = qs_kind%dudq_dftb3
568 53059330 : IF (PRESENT(elec_conf)) elec_conf => qs_kind%elec_conf
569 53059330 : IF (PRESENT(alpha_core_charge)) THEN
570 196413 : IF (ASSOCIATED(qs_kind%all_potential)) THEN
571 : CALL get_potential(potential=qs_kind%all_potential, &
572 44344 : alpha_core_charge=alpha_core_charge)
573 152069 : ELSE IF (ASSOCIATED(qs_kind%gth_potential)) THEN
574 : CALL get_potential(potential=qs_kind%gth_potential, &
575 150459 : alpha_core_charge=alpha_core_charge)
576 1610 : ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
577 : CALL get_potential(potential=qs_kind%sgp_potential, &
578 258 : alpha_core_charge=alpha_core_charge)
579 : ELSE
580 1352 : alpha_core_charge = 1.0_dp
581 : END IF
582 : END IF
583 53059330 : IF (PRESENT(ccore_charge)) THEN
584 80735 : IF (ASSOCIATED(qs_kind%all_potential)) THEN
585 : CALL get_potential(potential=qs_kind%all_potential, &
586 9106 : ccore_charge=ccore_charge)
587 71629 : ELSE IF (ASSOCIATED(qs_kind%gth_potential)) THEN
588 : CALL get_potential(potential=qs_kind%gth_potential, &
589 70721 : ccore_charge=ccore_charge)
590 908 : ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
591 : CALL get_potential(potential=qs_kind%sgp_potential, &
592 134 : ccore_charge=ccore_charge)
593 774 : ELSE IF (ASSOCIATED(qs_kind%upf_potential)) THEN
594 0 : CPABORT("UPF CCORE CHARGE RADIUS NOT AVAILABLE")
595 : ELSE
596 774 : ccore_charge = 0.0_dp
597 : END IF
598 : END IF
599 53059330 : IF (PRESENT(core_charge_radius)) THEN
600 79971 : IF (ASSOCIATED(qs_kind%all_potential)) THEN
601 : CALL get_potential(potential=qs_kind%all_potential, &
602 32078 : core_charge_radius=core_charge_radius)
603 47893 : ELSE IF (ASSOCIATED(qs_kind%gth_potential)) THEN
604 : CALL get_potential(potential=qs_kind%gth_potential, &
605 47449 : core_charge_radius=core_charge_radius)
606 444 : ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
607 : CALL get_potential(potential=qs_kind%sgp_potential, &
608 70 : core_charge_radius=core_charge_radius)
609 374 : ELSE IF (ASSOCIATED(qs_kind%upf_potential)) THEN
610 0 : CPABORT("UPF CORE CHARGE RADIUS NOT AVAILABLE")
611 : ELSE
612 374 : core_charge_radius = 0.0_dp
613 : END IF
614 : END IF
615 53059330 : IF (PRESENT(core_charge)) THEN
616 35170 : IF (ASSOCIATED(qs_kind%all_potential)) THEN
617 : CALL get_potential(potential=qs_kind%all_potential, &
618 365 : zeff=core_charge)
619 34805 : ELSE IF (ASSOCIATED(qs_kind%gth_potential)) THEN
620 : CALL get_potential(potential=qs_kind%gth_potential, &
621 34805 : zeff=core_charge)
622 0 : ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
623 : CALL get_potential(potential=qs_kind%sgp_potential, &
624 0 : zeff=core_charge)
625 0 : ELSE IF (ASSOCIATED(qs_kind%upf_potential)) THEN
626 0 : CPABORT("UPF CORE CHARGE NOT AVAILABLE")
627 : ELSE
628 0 : core_charge = 0.0_dp
629 : END IF
630 : END IF
631 :
632 53059330 : IF (PRESENT(zatom)) THEN
633 : ! Retrieve information on element
634 189230 : CALL get_ptable_info(qs_kind%element_symbol, ielement=zatom, found=found)
635 189230 : CPASSERT(found)
636 : END IF
637 :
638 53059330 : IF (PRESENT(zeff)) THEN
639 212757 : IF (ASSOCIATED(qs_kind%all_potential)) THEN
640 50961 : CALL get_potential(potential=qs_kind%all_potential, zeff=zeff)
641 161796 : ELSE IF (ASSOCIATED(qs_kind%gth_potential)) THEN
642 160578 : CALL get_potential(potential=qs_kind%gth_potential, zeff=zeff)
643 1218 : ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
644 254 : CALL get_potential(potential=qs_kind%sgp_potential, zeff=zeff)
645 964 : ELSE IF (ASSOCIATED(qs_kind%upf_potential)) THEN
646 65 : zeff = qs_kind%upf_potential%zion
647 : ELSE
648 899 : zeff = 0.0_dp
649 : END IF
650 : END IF
651 :
652 53059330 : IF (PRESENT(covalent_radius)) covalent_radius = qs_kind%covalent_radius
653 53059330 : IF (PRESENT(vdw_radius)) vdw_radius = qs_kind%vdw_radius
654 :
655 53059330 : IF (PRESENT(paw_proj_set)) paw_proj_set => qs_kind%paw_proj_set
656 53059330 : IF (PRESENT(paw_atom)) paw_atom = qs_kind%paw_atom
657 53059330 : IF (PRESENT(gpw_type_forced)) gpw_type_forced = qs_kind%gpw_type_forced
658 53059330 : IF (PRESENT(hard_radius)) hard_radius = qs_kind%hard_radius
659 53059330 : IF (PRESENT(hard0_radius)) hard0_radius = qs_kind%hard0_radius
660 53059330 : IF (PRESENT(max_rad_local)) max_rad_local = qs_kind%max_rad_local
661 53059330 : IF (PRESENT(harmonics)) harmonics => qs_kind%harmonics
662 53059330 : IF (PRESENT(max_s_harm)) THEN
663 7581383 : IF (ASSOCIATED(qs_kind%harmonics)) THEN
664 278488 : max_s_harm = qs_kind%harmonics%max_s_harm
665 : ELSE
666 7302895 : max_s_harm = 0
667 : END IF
668 : END IF
669 53059330 : IF (PRESENT(max_iso_not0)) THEN
670 7611757 : IF (ASSOCIATED(qs_kind%harmonics)) THEN
671 308862 : max_iso_not0 = qs_kind%harmonics%max_iso_not0
672 : ELSE
673 7302895 : max_iso_not0 = 0
674 : END IF
675 : END IF
676 53059330 : IF (PRESENT(grid_atom)) grid_atom => qs_kind%grid_atom
677 53059330 : IF (PRESENT(ngrid_ang)) ngrid_ang = qs_kind%ngrid_ang
678 53059330 : IF (PRESENT(ngrid_rad)) ngrid_rad = qs_kind%ngrid_rad
679 53059330 : IF (PRESENT(lmax_rho0)) lmax_rho0 = qs_kind%lmax_rho0
680 53059330 : IF (PRESENT(ghost)) ghost = qs_kind%ghost
681 53059330 : IF (PRESENT(floating)) floating = qs_kind%floating
682 53059330 : IF (PRESENT(dft_plus_u_atom)) dft_plus_u_atom = ASSOCIATED(qs_kind%dft_plus_u)
683 53059330 : IF (PRESENT(l_of_dft_plus_u)) THEN
684 4858 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
685 2418 : l_of_dft_plus_u = qs_kind%dft_plus_u%l
686 : ELSE
687 2440 : l_of_dft_plus_u = -1
688 : END IF
689 : END IF
690 53059330 : IF (PRESENT(n_of_dft_plus_u)) THEN
691 22 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
692 0 : n_of_dft_plus_u = qs_kind%dft_plus_u%n
693 : ELSE
694 22 : n_of_dft_plus_u = -1
695 : END IF
696 : END IF
697 53059330 : IF (PRESENT(u_minus_j)) THEN
698 4836 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
699 2418 : u_minus_j = qs_kind%dft_plus_u%u_minus_j
700 : ELSE
701 2418 : u_minus_j = 0.0_dp
702 : END IF
703 : END IF
704 53059330 : IF (PRESENT(u_minus_j_target)) THEN
705 4858 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
706 2418 : u_minus_j_target = qs_kind%dft_plus_u%u_minus_j_target
707 : ELSE
708 2440 : u_minus_j_target = 0.0_dp
709 : END IF
710 : END IF
711 53059330 : IF (PRESENT(U_of_dft_plus_u)) THEN
712 22 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
713 0 : U_of_dft_plus_u = qs_kind%dft_plus_u%U
714 : ELSE
715 22 : U_of_dft_plus_u = 0.0_dp
716 : END IF
717 : END IF
718 53059330 : IF (PRESENT(J_of_dft_plus_u)) THEN
719 22 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
720 0 : J_of_dft_plus_u = qs_kind%dft_plus_u%J
721 : ELSE
722 22 : J_of_dft_plus_u = 0.0_dp
723 : END IF
724 : END IF
725 53059330 : IF (PRESENT(alpha_of_dft_plus_u)) THEN
726 22 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
727 0 : alpha_of_dft_plus_u = qs_kind%dft_plus_u%alpha
728 : ELSE
729 22 : alpha_of_dft_plus_u = 0.0_dp
730 : END IF
731 : END IF
732 53059330 : IF (PRESENT(beta_of_dft_plus_u)) THEN
733 22 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
734 0 : beta_of_dft_plus_u = qs_kind%dft_plus_u%beta
735 : ELSE
736 22 : beta_of_dft_plus_u = 0.0_dp
737 : END IF
738 : END IF
739 53059330 : IF (PRESENT(J0_of_dft_plus_u)) THEN
740 22 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
741 0 : J0_of_dft_plus_u = qs_kind%dft_plus_u%J0
742 : ELSE
743 22 : J0_of_dft_plus_u = 0.0_dp
744 : END IF
745 : END IF
746 53059330 : IF (PRESENT(occupation_of_dft_plus_u)) THEN
747 22 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
748 0 : occupation_of_dft_plus_u = qs_kind%dft_plus_u%occupation
749 : ELSE
750 22 : occupation_of_dft_plus_u = -1.0_dp
751 : END IF
752 : END IF
753 :
754 53059330 : IF (PRESENT(init_u_ramping_each_scf)) THEN
755 160 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
756 80 : init_u_ramping_each_scf = qs_kind%dft_plus_u%init_u_ramping_each_scf
757 : ELSE
758 80 : init_u_ramping_each_scf = .FALSE.
759 : END IF
760 : END IF
761 53059330 : IF (PRESENT(u_ramping)) THEN
762 4996 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
763 2498 : u_ramping = qs_kind%dft_plus_u%u_ramping
764 : ELSE
765 2498 : u_ramping = 0.0_dp
766 : END IF
767 : END IF
768 53059330 : IF (PRESENT(eps_u_ramping)) THEN
769 4836 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
770 2418 : eps_u_ramping = qs_kind%dft_plus_u%eps_u_ramping
771 : ELSE
772 2418 : eps_u_ramping = 1.0E-5_dp
773 : END IF
774 : END IF
775 53059330 : IF (PRESENT(nelec)) THEN
776 3640 : NULLIFY (nelec)
777 3640 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
778 1820 : IF (ASSOCIATED(qs_kind%dft_plus_u%nelec)) THEN
779 0 : nelec => qs_kind%dft_plus_u%nelec
780 : END IF
781 : END IF
782 : END IF
783 53059330 : IF (PRESENT(orbitals)) THEN
784 3912 : NULLIFY (orbitals)
785 3912 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
786 1956 : IF (ASSOCIATED(qs_kind%dft_plus_u%orbitals)) THEN
787 112 : orbitals => qs_kind%dft_plus_u%orbitals
788 : END IF
789 : END IF
790 : END IF
791 53059330 : IF (PRESENT(eps_scf)) THEN
792 3912 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
793 1956 : eps_scf = qs_kind%dft_plus_u%eps_scf
794 : ELSE
795 1956 : eps_scf = 1.0E30_dp
796 : END IF
797 : END IF
798 53059330 : IF (PRESENT(max_scf)) THEN
799 3912 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
800 1956 : max_scf = qs_kind%dft_plus_u%max_scf
801 : ELSE
802 1956 : max_scf = -1
803 : END IF
804 : END IF
805 53059330 : IF (PRESENT(smear)) THEN
806 3912 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
807 1956 : smear = qs_kind%dft_plus_u%smear
808 : ELSE
809 1956 : smear = .FALSE.
810 : END IF
811 : END IF
812 53059330 : IF (PRESENT(dispersion)) dispersion => qs_kind%dispersion
813 53059330 : IF (PRESENT(bs_occupation)) bs_occupation = qs_kind%bs_occupation
814 53059330 : IF (PRESENT(addel)) addel => qs_kind%addel
815 53059330 : IF (PRESENT(laddel)) laddel => qs_kind%laddel
816 53059330 : IF (PRESENT(naddel)) naddel => qs_kind%naddel
817 :
818 53059330 : IF (PRESENT(magnetization)) magnetization = qs_kind%magnetization
819 :
820 53059330 : IF (PRESENT(no_optimize)) no_optimize = qs_kind%no_optimize
821 :
822 53059330 : IF (PRESENT(reltmat)) reltmat => qs_kind%reltmat
823 :
824 53059330 : IF (PRESENT(mao)) mao = qs_kind%mao
825 :
826 53059330 : IF (PRESENT(lmax_dftb)) lmax_dftb = qs_kind%lmax_dftb
827 :
828 53059330 : IF (PRESENT(pao_basis_size)) pao_basis_size = qs_kind%pao_basis_size
829 53059330 : IF (PRESENT(pao_model_file)) pao_model_file = qs_kind%pao_model_file
830 53059330 : IF (PRESENT(pao_potentials)) pao_potentials => qs_kind%pao_potentials
831 53059330 : IF (PRESENT(pao_descriptors)) pao_descriptors => qs_kind%pao_descriptors
832 53059330 : END SUBROUTINE get_qs_kind
833 :
834 : ! **************************************************************************************************
835 : !> \brief Get attributes of an atomic kind set.
836 : !> \param qs_kind_set ...
837 : !> \param all_potential_present ...
838 : !> \param tnadd_potential_present ...
839 : !> \param gth_potential_present ...
840 : !> \param sgp_potential_present ...
841 : !> \param paw_atom_present ...
842 : !> \param dft_plus_u_atom_present ...
843 : !> \param maxcgf ...
844 : !> \param maxsgf ...
845 : !> \param maxco ...
846 : !> \param maxco_proj ...
847 : !> \param maxgtops ...
848 : !> \param maxlgto ...
849 : !> \param maxlprj ...
850 : !> \param maxnset ...
851 : !> \param maxsgf_set ...
852 : !> \param ncgf ...
853 : !> \param npgf ...
854 : !> \param nset ...
855 : !> \param nsgf ...
856 : !> \param nshell ...
857 : !> \param maxpol ...
858 : !> \param maxlppl ...
859 : !> \param maxlppnl ...
860 : !> \param maxppnl ...
861 : !> \param nelectron ...
862 : !> \param maxder ...
863 : !> \param max_ngrid_rad ...
864 : !> \param max_sph_harm ...
865 : !> \param maxg_iso_not0 ...
866 : !> \param lmax_rho0 ...
867 : !> \param basis_rcut ...
868 : !> \param basis_type ...
869 : !> \param total_zeff_corr ... [SGh]
870 : ! **************************************************************************************************
871 3571364 : SUBROUTINE get_qs_kind_set(qs_kind_set, &
872 : all_potential_present, tnadd_potential_present, gth_potential_present, &
873 : sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, &
874 : maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, &
875 : ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, &
876 : nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, &
877 : basis_rcut, &
878 : basis_type, total_zeff_corr)
879 :
880 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
881 : LOGICAL, INTENT(OUT), OPTIONAL :: all_potential_present, tnadd_potential_present, &
882 : gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present
883 : INTEGER, INTENT(OUT), OPTIONAL :: maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, &
884 : maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, &
885 : maxppnl, nelectron
886 : INTEGER, INTENT(IN), OPTIONAL :: maxder
887 : INTEGER, INTENT(OUT), OPTIONAL :: max_ngrid_rad, max_sph_harm, &
888 : maxg_iso_not0, lmax_rho0
889 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: basis_rcut
890 : CHARACTER(len=*), OPTIONAL :: basis_type
891 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: total_zeff_corr
892 :
893 : CHARACTER(len=default_string_length) :: my_basis_type
894 : INTEGER :: ikind, imax, lmax_rho0_kind, &
895 : max_iso_not0, max_s_harm, n, &
896 : ngrid_rad, nkind, nrloc(10), &
897 : nrpot(1:15, 0:10)
898 : LOGICAL :: dft_plus_u_atom, ecp_semi_local, paw_atom
899 : REAL(KIND=dp) :: brcut, zeff, zeff_correction
900 : TYPE(all_potential_type), POINTER :: all_potential
901 : TYPE(gth_potential_type), POINTER :: gth_potential
902 : TYPE(gto_basis_set_type), POINTER :: tmp_basis_set
903 : TYPE(local_potential_type), POINTER :: tnadd_potential
904 : TYPE(paw_proj_set_type), POINTER :: paw_proj_set
905 : TYPE(qs_dftb_atom_type), POINTER :: dftb_parameter
906 : TYPE(qs_kind_type), POINTER :: qs_kind
907 : TYPE(sgp_potential_type), POINTER :: sgp_potential
908 :
909 3571364 : IF (PRESENT(basis_type)) THEN
910 3317882 : my_basis_type = basis_type
911 : ELSE
912 253482 : my_basis_type = "ORB"
913 : END IF
914 :
915 3571364 : IF (ASSOCIATED(qs_kind_set)) THEN
916 :
917 3571364 : IF (PRESENT(maxcgf)) maxcgf = 0
918 3571364 : IF (PRESENT(maxco)) maxco = 0
919 3571364 : IF (PRESENT(maxco_proj)) maxco_proj = 0
920 3571364 : IF (PRESENT(maxg_iso_not0)) maxg_iso_not0 = 0
921 3571364 : IF (PRESENT(maxgtops)) maxgtops = 0
922 3571364 : IF (PRESENT(maxlgto)) maxlgto = -1
923 3571364 : IF (PRESENT(maxlppl)) maxlppl = -1
924 3571364 : IF (PRESENT(maxlppnl)) maxlppnl = -1
925 3571364 : IF (PRESENT(maxpol)) maxpol = -1
926 3571364 : IF (PRESENT(maxlprj)) maxlprj = -1
927 3571364 : IF (PRESENT(maxnset)) maxnset = 0
928 3571364 : IF (PRESENT(maxppnl)) maxppnl = 0
929 3571364 : IF (PRESENT(maxsgf)) maxsgf = 0
930 3571364 : IF (PRESENT(maxsgf_set)) maxsgf_set = 0
931 3571364 : IF (PRESENT(ncgf)) ncgf = 0
932 3571364 : IF (PRESENT(nelectron)) nelectron = 0
933 3571364 : IF (PRESENT(npgf)) npgf = 0
934 3571364 : IF (PRESENT(nset)) nset = 0
935 3571364 : IF (PRESENT(nsgf)) nsgf = 0
936 3571364 : IF (PRESENT(nshell)) nshell = 0
937 3571364 : IF (PRESENT(all_potential_present)) all_potential_present = .FALSE.
938 3571364 : IF (PRESENT(tnadd_potential_present)) tnadd_potential_present = .FALSE.
939 3571364 : IF (PRESENT(gth_potential_present)) gth_potential_present = .FALSE.
940 3571364 : IF (PRESENT(sgp_potential_present)) sgp_potential_present = .FALSE.
941 3571364 : IF (PRESENT(paw_atom_present)) paw_atom_present = .FALSE.
942 3571364 : IF (PRESENT(max_ngrid_rad)) max_ngrid_rad = 0
943 3571364 : IF (PRESENT(max_sph_harm)) max_sph_harm = 0
944 3571364 : IF (PRESENT(lmax_rho0)) lmax_rho0 = 0
945 3571364 : IF (PRESENT(basis_rcut)) basis_rcut = 0.0_dp
946 3571364 : IF (PRESENT(total_zeff_corr)) total_zeff_corr = 0.0_dp
947 :
948 3571364 : nkind = SIZE(qs_kind_set)
949 11152747 : DO ikind = 1, nkind
950 7581383 : qs_kind => qs_kind_set(ikind)
951 : CALL get_qs_kind(qs_kind=qs_kind, &
952 : all_potential=all_potential, &
953 : tnadd_potential=tnadd_potential, &
954 : gth_potential=gth_potential, &
955 : sgp_potential=sgp_potential, &
956 : paw_proj_set=paw_proj_set, &
957 : dftb_parameter=dftb_parameter, &
958 : ngrid_rad=ngrid_rad, &
959 : max_s_harm=max_s_harm, &
960 : max_iso_not0=max_iso_not0, &
961 : paw_atom=paw_atom, &
962 : dft_plus_u_atom=dft_plus_u_atom, &
963 7581383 : lmax_rho0=lmax_rho0_kind)
964 :
965 7581383 : IF (PRESENT(maxlppl) .AND. ASSOCIATED(gth_potential)) THEN
966 42874 : CALL get_potential(potential=gth_potential, nexp_ppl=n)
967 42874 : maxlppl = MAX(maxlppl, 2*(n - 1))
968 7185 : ELSEIF (PRESENT(maxlppl) .AND. ASSOCIATED(sgp_potential)) THEN
969 74 : CALL get_potential(potential=sgp_potential, nrloc=nrloc, ecp_semi_local=ecp_semi_local)
970 814 : n = MAXVAL(nrloc) - 2
971 74 : maxlppl = MAX(maxlppl, 2*(n - 1))
972 74 : IF (ecp_semi_local) THEN
973 44 : CALL get_potential(potential=sgp_potential, sl_lmax=imax, nrpot=nrpot)
974 7788 : n = MAXVAL(nrpot) - 2
975 44 : n = 2*(n - 1) + imax
976 44 : maxlppl = MAX(maxlppl, n)
977 : END IF
978 : END IF
979 :
980 7581383 : IF (PRESENT(maxlppnl) .AND. ASSOCIATED(gth_potential)) THEN
981 39885 : CALL get_potential(potential=gth_potential, lprj_ppnl_max=imax)
982 39885 : maxlppnl = MAX(maxlppnl, imax)
983 7147 : ELSEIF (PRESENT(maxlppnl) .AND. ASSOCIATED(sgp_potential)) THEN
984 46 : CALL get_potential(potential=sgp_potential, lmax=imax)
985 46 : maxlppnl = MAX(maxlppnl, imax)
986 : END IF
987 :
988 7581383 : IF (PRESENT(maxpol) .AND. ASSOCIATED(tnadd_potential)) THEN
989 66 : CALL get_potential(potential=tnadd_potential, npol=n)
990 66 : maxpol = MAX(maxpol, 2*(n - 1))
991 : END IF
992 :
993 7581383 : IF (PRESENT(maxco_proj) .AND. ASSOCIATED(paw_proj_set)) THEN
994 4204 : CALL get_paw_proj_set(paw_proj_set=paw_proj_set, ncgauprj=imax)
995 4204 : maxco_proj = MAX(maxco_proj, imax)
996 : END IF
997 :
998 7581383 : IF (PRESENT(maxlprj) .AND. ASSOCIATED(paw_proj_set)) THEN
999 4204 : CALL get_paw_proj_set(paw_proj_set=paw_proj_set, maxl=imax)
1000 4204 : maxlprj = MAX(maxlprj, imax)
1001 : END IF
1002 :
1003 7581383 : IF (PRESENT(maxppnl) .AND. ASSOCIATED(gth_potential)) THEN
1004 27516 : CALL get_potential(potential=gth_potential, nppnl=imax)
1005 27516 : maxppnl = MAX(maxppnl, imax)
1006 232 : ELSEIF (PRESENT(maxppnl) .AND. ASSOCIATED(sgp_potential)) THEN
1007 10 : CALL get_potential(potential=sgp_potential, nppnl=imax)
1008 10 : maxppnl = MAX(maxppnl, imax)
1009 : END IF
1010 :
1011 : CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis_set, &
1012 7581383 : basis_type=my_basis_type)
1013 :
1014 7581383 : IF (PRESENT(maxcgf)) THEN
1015 0 : IF (ASSOCIATED(tmp_basis_set)) THEN
1016 0 : CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, ncgf=imax)
1017 0 : maxcgf = MAX(maxcgf, imax)
1018 0 : ELSE IF (ASSOCIATED(qs_kind%dftb_parameter)) THEN
1019 0 : CALL get_dftb_atom_param(dftb_parameter=dftb_parameter, lmax=imax)
1020 0 : imax = ((imax + 1)*(imax + 2)*(imax + 3))/6
1021 0 : maxcgf = MAX(maxcgf, imax)
1022 : END IF
1023 : END IF
1024 :
1025 7581383 : IF (PRESENT(maxco)) THEN
1026 6949329 : IF (ASSOCIATED(tmp_basis_set)) THEN
1027 6949321 : IF (PRESENT(maxder)) THEN
1028 : CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, &
1029 0 : maxco=imax, maxder=maxder)
1030 : ELSE
1031 6949321 : CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, maxco=imax)
1032 : END IF
1033 6949321 : maxco = MAX(maxco, imax)
1034 : END IF
1035 6949329 : IF (ASSOCIATED(gth_potential)) THEN
1036 630933 : CALL get_potential(potential=gth_potential, lprj_ppnl_max=imax)
1037 630933 : maxco = MAX(maxco, ncoset(imax))
1038 : END IF
1039 6949329 : IF (ASSOCIATED(sgp_potential)) THEN
1040 584 : CALL get_potential(potential=sgp_potential, lmax=imax)
1041 584 : maxco = MAX(maxco, ncoset(imax))
1042 584 : CALL get_potential(potential=sgp_potential, sl_lmax=imax)
1043 584 : maxco = MAX(maxco, ncoset(imax))
1044 : END IF
1045 : END IF
1046 :
1047 7581383 : IF (PRESENT(maxgtops)) THEN
1048 95390 : IF (ASSOCIATED(tmp_basis_set)) THEN
1049 95390 : CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, maxso=imax, nset=n)
1050 95390 : maxgtops = MAX(maxgtops, n*imax)
1051 : END IF
1052 : END IF
1053 :
1054 7581383 : IF (PRESENT(maxlgto)) THEN
1055 6552480 : IF (ASSOCIATED(tmp_basis_set)) THEN
1056 6531873 : CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, maxl=imax)
1057 6531873 : maxlgto = MAX(maxlgto, imax)
1058 20607 : ELSE IF (ASSOCIATED(qs_kind%dftb_parameter)) THEN
1059 2516 : CALL get_dftb_atom_param(dftb_parameter=dftb_parameter, lmax=imax)
1060 2516 : maxlgto = MAX(maxlgto, imax)
1061 : END IF
1062 : END IF
1063 :
1064 7581383 : IF (PRESENT(maxnset)) THEN
1065 73927 : IF (ASSOCIATED(tmp_basis_set)) THEN
1066 73927 : CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, nset=n)
1067 73927 : maxnset = MAX(maxnset, n)
1068 : END IF
1069 : END IF
1070 :
1071 7581383 : IF (PRESENT(maxsgf)) THEN
1072 6659320 : IF (ASSOCIATED(tmp_basis_set)) THEN
1073 6659296 : CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, nsgf=imax)
1074 6659296 : maxsgf = MAX(maxsgf, imax)
1075 : END IF
1076 : END IF
1077 :
1078 7581383 : IF (PRESENT(maxsgf_set)) THEN
1079 437329 : IF (ASSOCIATED(tmp_basis_set)) THEN
1080 437329 : CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, maxsgf_set=imax)
1081 437329 : maxsgf_set = MAX(maxsgf_set, imax)
1082 : END IF
1083 : END IF
1084 :
1085 7581383 : IF (PRESENT(ncgf)) THEN
1086 31998 : IF (ASSOCIATED(tmp_basis_set)) THEN
1087 12934 : CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, ncgf=n)
1088 12934 : ncgf = ncgf + n*qs_kind_set(ikind)%natom
1089 19064 : ELSE IF (ASSOCIATED(qs_kind%dftb_parameter)) THEN
1090 987 : CALL get_dftb_atom_param(dftb_parameter=dftb_parameter, lmax=imax)
1091 987 : n = ((imax + 1)*(imax + 2)*(imax + 3))/6
1092 987 : ncgf = ncgf + n*qs_kind_set(ikind)%natom
1093 : END IF
1094 : END IF
1095 :
1096 7581383 : IF (PRESENT(npgf)) THEN
1097 25660 : IF (ASSOCIATED(tmp_basis_set)) THEN
1098 6623 : CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, npgf_sum=n)
1099 6623 : npgf = npgf + n*qs_kind_set(ikind)%natom
1100 : END IF
1101 : END IF
1102 :
1103 7581383 : IF (PRESENT(nset)) THEN
1104 25660 : IF (ASSOCIATED(tmp_basis_set)) THEN
1105 6623 : CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, nset=n)
1106 6623 : nset = nset + n*qs_kind_set(ikind)%natom
1107 : END IF
1108 : END IF
1109 :
1110 7581383 : IF (PRESENT(nsgf)) THEN
1111 95030 : IF (ASSOCIATED(tmp_basis_set)) THEN
1112 62096 : CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, nsgf=n)
1113 62096 : nsgf = nsgf + n*qs_kind_set(ikind)%natom
1114 32934 : ELSE IF (ASSOCIATED(qs_kind%dftb_parameter)) THEN
1115 14855 : CALL get_dftb_atom_param(dftb_parameter=dftb_parameter, natorb=n)
1116 14855 : nsgf = nsgf + n*qs_kind_set(ikind)%natom
1117 : END IF
1118 : END IF
1119 :
1120 7581383 : IF (PRESENT(nshell)) THEN
1121 25660 : IF (ASSOCIATED(tmp_basis_set)) THEN
1122 6623 : CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, nshell_sum=n)
1123 6623 : nshell = nshell + n*qs_kind_set(ikind)%natom
1124 19037 : ELSE IF (ASSOCIATED(qs_kind%dftb_parameter)) THEN
1125 960 : CALL get_dftb_atom_param(dftb_parameter=dftb_parameter, lmax=n)
1126 960 : nshell = nshell + (n + 1)*qs_kind_set(ikind)%natom
1127 : END IF
1128 : END IF
1129 :
1130 7581383 : IF (PRESENT(nelectron)) THEN
1131 197542 : IF (ASSOCIATED(qs_kind%all_potential)) THEN
1132 : CALL get_potential(potential=qs_kind%all_potential, &
1133 16670 : zeff=zeff, zeff_correction=zeff_correction)
1134 180872 : ELSE IF (ASSOCIATED(qs_kind%gth_potential)) THEN
1135 : CALL get_potential(potential=qs_kind%gth_potential, &
1136 179690 : zeff=zeff, zeff_correction=zeff_correction)
1137 1182 : ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
1138 : CALL get_potential(potential=qs_kind%sgp_potential, &
1139 396 : zeff=zeff, zeff_correction=zeff_correction)
1140 : ELSE
1141 786 : zeff = 0.0_dp
1142 786 : zeff_correction = 0.0_dp
1143 : END IF
1144 197542 : nelectron = nelectron + qs_kind_set(ikind)%natom*NINT(zeff - zeff_correction)
1145 : END IF
1146 :
1147 7581383 : IF (PRESENT(basis_rcut)) THEN
1148 234 : IF (ASSOCIATED(tmp_basis_set)) THEN
1149 0 : CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, kind_radius=brcut)
1150 0 : basis_rcut = MAX(basis_rcut, brcut)
1151 234 : ELSE IF (ASSOCIATED(qs_kind%dftb_parameter)) THEN
1152 234 : CALL get_dftb_atom_param(dftb_parameter=dftb_parameter, cutoff=brcut)
1153 234 : basis_rcut = MAX(basis_rcut, brcut)
1154 : END IF
1155 : END IF
1156 :
1157 7581383 : IF (PRESENT(total_zeff_corr)) THEN
1158 12745 : IF (ASSOCIATED(qs_kind%all_potential)) THEN
1159 : CALL get_potential(potential=qs_kind%all_potential, &
1160 4434 : zeff=zeff, zeff_correction=zeff_correction)
1161 8311 : ELSE IF (ASSOCIATED(qs_kind%gth_potential)) THEN
1162 : CALL get_potential(potential=qs_kind%gth_potential, &
1163 8143 : zeff=zeff, zeff_correction=zeff_correction)
1164 168 : ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
1165 : CALL get_potential(potential=qs_kind%sgp_potential, &
1166 24 : zeff=zeff, zeff_correction=zeff_correction)
1167 : ELSE
1168 144 : zeff = 0.0_dp
1169 144 : zeff_correction = 0.0_dp
1170 : END IF
1171 12745 : total_zeff_corr = total_zeff_corr + qs_kind_set(ikind)%natom*zeff_correction
1172 : END IF
1173 :
1174 7581383 : IF (PRESENT(all_potential_present)) THEN
1175 60513 : IF (ASSOCIATED(all_potential)) THEN
1176 34808 : all_potential_present = .TRUE.
1177 : END IF
1178 : END IF
1179 :
1180 7581383 : IF (PRESENT(tnadd_potential_present)) THEN
1181 0 : IF (ASSOCIATED(tnadd_potential)) THEN
1182 0 : tnadd_potential_present = .TRUE.
1183 : END IF
1184 : END IF
1185 :
1186 7581383 : IF (PRESENT(gth_potential_present)) THEN
1187 48386 : IF (ASSOCIATED(gth_potential)) THEN
1188 17852 : gth_potential_present = .TRUE.
1189 : END IF
1190 : END IF
1191 :
1192 7581383 : IF (PRESENT(sgp_potential_present)) THEN
1193 48386 : IF (ASSOCIATED(sgp_potential)) THEN
1194 32 : sgp_potential_present = .TRUE.
1195 : END IF
1196 : END IF
1197 :
1198 7581383 : IF (PRESENT(paw_atom_present)) THEN
1199 47716 : IF (paw_atom) THEN
1200 2914 : paw_atom_present = .TRUE.
1201 : END IF
1202 : END IF
1203 :
1204 7581383 : IF (PRESENT(dft_plus_u_atom_present)) THEN
1205 12745 : IF (dft_plus_u_atom) THEN
1206 32 : dft_plus_u_atom_present = .TRUE.
1207 : END IF
1208 : END IF
1209 :
1210 7581383 : IF (PRESENT(max_ngrid_rad)) THEN
1211 0 : max_ngrid_rad = MAX(max_ngrid_rad, ngrid_rad)
1212 : END IF
1213 :
1214 7581383 : IF (PRESENT(max_sph_harm)) THEN
1215 0 : max_sph_harm = MAX(max_sph_harm, max_s_harm)
1216 : END IF
1217 :
1218 7581383 : IF (PRESENT(maxg_iso_not0)) THEN
1219 30374 : maxg_iso_not0 = MAX(maxg_iso_not0, max_iso_not0)
1220 : END IF
1221 :
1222 18734130 : IF (PRESENT(lmax_rho0)) THEN
1223 0 : lmax_rho0 = MAX(lmax_rho0, lmax_rho0_kind)
1224 : END IF
1225 :
1226 : END DO
1227 : ELSE
1228 0 : CPABORT("The pointer qs_kind_set is not associated")
1229 : END IF
1230 :
1231 3571364 : END SUBROUTINE get_qs_kind_set
1232 :
1233 : ! **************************************************************************************************
1234 : !> \brief Initialise an atomic kind data set.
1235 : !> \param qs_kind ...
1236 : !> \author Creation (11.01.2002,MK)
1237 : !> 20.09.2002 adapted for pol/kg use, gtb
1238 : ! **************************************************************************************************
1239 12745 : SUBROUTINE init_qs_kind(qs_kind)
1240 : TYPE(qs_kind_type), POINTER :: qs_kind
1241 :
1242 : CHARACTER(len=*), PARAMETER :: routineN = 'init_qs_kind'
1243 :
1244 : CHARACTER(LEN=default_string_length) :: basis_type
1245 : INTEGER :: handle, i
1246 : TYPE(gto_basis_set_type), POINTER :: tmp_basis_set
1247 :
1248 12745 : CALL timeset(routineN, handle)
1249 :
1250 12745 : CPASSERT(ASSOCIATED(qs_kind))
1251 :
1252 12745 : IF (ASSOCIATED(qs_kind%gth_potential)) THEN
1253 8143 : CALL init_potential(qs_kind%gth_potential)
1254 4602 : ELSEIF (ASSOCIATED(qs_kind%sgp_potential)) THEN
1255 24 : CALL init_potential(qs_kind%sgp_potential)
1256 : END IF
1257 :
1258 267645 : DO i = 1, SIZE(qs_kind%basis_sets, 1)
1259 254900 : NULLIFY (tmp_basis_set)
1260 : CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis_set, &
1261 254900 : inumbas=i, basis_type=basis_type)
1262 254900 : IF (basis_type == "") CYCLE
1263 27351 : IF (basis_type == "AUX") THEN
1264 0 : IF (tmp_basis_set%norm_type < 0) tmp_basis_set%norm_type = 1
1265 0 : CALL init_aux_basis_set(tmp_basis_set)
1266 : ELSE
1267 14606 : IF (tmp_basis_set%norm_type < 0) tmp_basis_set%norm_type = 2
1268 14606 : CALL init_orb_basis_set(tmp_basis_set)
1269 : END IF
1270 : END DO
1271 :
1272 12745 : CALL timestop(handle)
1273 :
1274 12745 : END SUBROUTINE init_qs_kind
1275 :
1276 : ! **************************************************************************************************
1277 : !> \brief Initialise an atomic kind set data set.
1278 : !> \param qs_kind_set ...
1279 : !> \author - Creation (17.01.2002,MK)
1280 : !> - 20.09.2002 para_env passed (gt)
1281 : ! **************************************************************************************************
1282 6686 : SUBROUTINE init_qs_kind_set(qs_kind_set)
1283 :
1284 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1285 :
1286 : CHARACTER(len=*), PARAMETER :: routineN = 'init_qs_kind_set'
1287 :
1288 : INTEGER :: handle, ikind
1289 : TYPE(qs_kind_type), POINTER :: qs_kind
1290 :
1291 6686 : CALL timeset(routineN, handle)
1292 :
1293 6686 : IF (.NOT. ASSOCIATED(qs_kind_set)) THEN
1294 0 : CPABORT("init_qs_kind_set: The pointer qs_kind_set is not associated")
1295 : END IF
1296 :
1297 19431 : DO ikind = 1, SIZE(qs_kind_set)
1298 12745 : qs_kind => qs_kind_set(ikind)
1299 19431 : CALL init_qs_kind(qs_kind)
1300 : END DO
1301 :
1302 6686 : CALL timestop(handle)
1303 :
1304 6686 : END SUBROUTINE init_qs_kind_set
1305 :
1306 : ! **************************************************************************************************
1307 : !> \brief ...
1308 : !> \param qs_kind_set ...
1309 : !> \param qs_control ...
1310 : !> \param force_env_section ...
1311 : !> \param modify_qs_control whether the qs_control should be modified
1312 : ! **************************************************************************************************
1313 1004 : SUBROUTINE init_gapw_basis_set(qs_kind_set, qs_control, force_env_section, modify_qs_control)
1314 :
1315 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1316 : TYPE(qs_control_type), POINTER :: qs_control
1317 : TYPE(section_vals_type), POINTER :: force_env_section
1318 : LOGICAL, OPTIONAL :: modify_qs_control
1319 :
1320 : CHARACTER(LEN=default_string_length) :: bsname
1321 : INTEGER :: bas1c, ikind, ilevel, nkind
1322 : LOGICAL :: gpw, my_mod_control, paw_atom
1323 : REAL(dp) :: max_rad_local_type, rc
1324 : TYPE(gto_basis_set_type), POINTER :: basis_1c, orb_basis, soft_basis
1325 : TYPE(paw_proj_set_type), POINTER :: paw_proj
1326 : TYPE(qs_kind_type), POINTER :: qs_kind
1327 :
1328 1004 : my_mod_control = .TRUE.
1329 1004 : IF (PRESENT(modify_qs_control)) THEN
1330 84 : my_mod_control = modify_qs_control
1331 : END IF
1332 :
1333 1004 : IF (ASSOCIATED(qs_kind_set)) THEN
1334 :
1335 1004 : IF (my_mod_control) qs_control%gapw_control%non_paw_atoms = .FALSE.
1336 1004 : nkind = SIZE(qs_kind_set)
1337 :
1338 2932 : DO ikind = 1, nkind
1339 :
1340 1928 : qs_kind => qs_kind_set(ikind)
1341 :
1342 1928 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis)
1343 : CALL get_qs_kind(qs_kind=qs_kind, hard_radius=rc, &
1344 1928 : max_rad_local=max_rad_local_type, gpw_type_forced=gpw)
1345 :
1346 1928 : NULLIFY (soft_basis)
1347 1928 : CALL allocate_gto_basis_set(soft_basis)
1348 : CALL create_soft_basis(orb_basis, soft_basis, &
1349 : qs_control%gapw_control%eps_fit, rc, paw_atom, &
1350 1928 : qs_control%gapw_control%force_paw, gpw)
1351 1928 : CALL add_basis_set_to_container(qs_kind%basis_sets, soft_basis, "ORB_SOFT")
1352 1928 : CALL set_qs_kind(qs_kind=qs_kind, paw_atom=paw_atom)
1353 :
1354 1928 : bas1c = qs_control%gapw_control%basis_1c
1355 1928 : NULLIFY (basis_1c)
1356 1886 : SELECT CASE (bas1c)
1357 : CASE (gapw_1c_orb)
1358 1886 : ilevel = 0
1359 : CASE (gapw_1c_small)
1360 26 : ilevel = 1
1361 : CASE (gapw_1c_medium)
1362 4 : ilevel = 2
1363 : CASE (gapw_1c_large)
1364 8 : ilevel = 3
1365 : CASE (gapw_1c_very_large)
1366 4 : ilevel = 4
1367 : CASE DEFAULT
1368 1928 : CPABORT("basis_1c type")
1369 : END SELECT
1370 1928 : CALL remove_basis_from_container(qs_kind%basis_sets, basis_type="GAPW_1C")
1371 1928 : CALL create_1c_basis(orb_basis, soft_basis, basis_1c, ilevel)
1372 1928 : CALL get_gto_basis_set(gto_basis_set=orb_basis, name=bsname)
1373 1928 : basis_1c%name = TRIM(bsname)//"_1c"
1374 1928 : CALL add_basis_set_to_container(qs_kind%basis_sets, basis_1c, "GAPW_1C")
1375 1928 : IF (paw_atom) THEN
1376 1618 : CALL allocate_paw_proj_set(qs_kind%paw_proj_set)
1377 1618 : CALL get_qs_kind(qs_kind=qs_kind, paw_proj_set=paw_proj)
1378 : CALL projectors(paw_proj, basis_1c, orb_basis, rc, qs_control, &
1379 1618 : max_rad_local_type, force_env_section)
1380 : ELSE
1381 310 : IF (my_mod_control) qs_control%gapw_control%non_paw_atoms = .TRUE.
1382 : END IF
1383 :
1384 : ! grid_atom and harmonics are allocated even if NOT PAW_ATOM
1385 1928 : NULLIFY (qs_kind%grid_atom, qs_kind%harmonics)
1386 1928 : CALL allocate_grid_atom(qs_kind%grid_atom)
1387 6788 : CALL allocate_harmonics_atom(qs_kind%harmonics)
1388 :
1389 : END DO
1390 :
1391 1004 : IF (my_mod_control) THEN
1392 920 : IF (qs_control%gapw_control%non_paw_atoms) THEN
1393 150 : qs_control%gapw_control%nopaw_as_gpw = .TRUE.
1394 : ELSE
1395 770 : qs_control%gapw_control%nopaw_as_gpw = .FALSE.
1396 : END IF
1397 : END IF
1398 : ELSE
1399 0 : CPABORT("The pointer qs_kind_set is not associated")
1400 : END IF
1401 :
1402 1004 : END SUBROUTINE init_gapw_basis_set
1403 : ! **************************************************************************************************
1404 : !> \brief ...
1405 : !> \param qs_kind_set ...
1406 : ! **************************************************************************************************
1407 1004 : SUBROUTINE init_gapw_nlcc(qs_kind_set)
1408 :
1409 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1410 :
1411 : INTEGER :: i, ic, ikind, n_nlcc, nc, nexp_nlcc, &
1412 : nkind, nr
1413 1004 : INTEGER, DIMENSION(:), POINTER :: nct_nlcc
1414 : LOGICAL :: nlcc, nlcc_type, paw_atom
1415 : REAL(dp) :: alpha, coa, cval
1416 1004 : REAL(KIND=dp), DIMENSION(:), POINTER :: a_nlcc, alpha_nlcc, c_nlcc, fe, rc, rr
1417 1004 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cval_nlcc, den
1418 : TYPE(gth_potential_type), POINTER :: gth_potential
1419 : TYPE(qs_kind_type), POINTER :: qs_kind
1420 : TYPE(sgp_potential_type), POINTER :: sgp_potential
1421 :
1422 1004 : IF (ASSOCIATED(qs_kind_set)) THEN
1423 1004 : nlcc = has_nlcc(qs_kind_set)
1424 1004 : IF (nlcc) THEN
1425 2 : nkind = SIZE(qs_kind_set)
1426 4 : DO ikind = 1, nkind
1427 2 : qs_kind => qs_kind_set(ikind)
1428 2 : CALL get_qs_kind(qs_kind, paw_atom=paw_atom)
1429 4 : IF (paw_atom) THEN
1430 2 : CALL get_qs_kind(qs_kind, gth_potential=gth_potential)
1431 2 : CALL get_qs_kind(qs_kind, sgp_potential=sgp_potential)
1432 2 : IF (ASSOCIATED(gth_potential)) THEN
1433 : CALL get_potential(potential=gth_potential, nlcc_present=nlcc_type, &
1434 2 : nexp_nlcc=nexp_nlcc, alpha_nlcc=alpha_nlcc, nct_nlcc=nct_nlcc, cval_nlcc=cval_nlcc)
1435 2 : IF (nlcc_type) THEN
1436 2 : nr = qs_kind%grid_atom%nr
1437 2 : rr => qs_kind%grid_atom%rad
1438 12 : ALLOCATE (qs_kind%nlcc_pot(nr, 2), rc(nr), fe(nr))
1439 6 : den => qs_kind%nlcc_pot
1440 206 : den = 0.0_dp
1441 4 : DO i = 1, nexp_nlcc
1442 2 : alpha = alpha_nlcc(i)
1443 202 : rc(:) = rr(:)/alpha
1444 202 : fe(:) = EXP(-0.5_dp*rc(:)*rc(:))
1445 2 : nc = nct_nlcc(i)
1446 8 : DO ic = 1, nc
1447 4 : cval = cval_nlcc(ic, i)
1448 4 : coa = cval/alpha
1449 404 : den(:, 1) = den(:, 1) + fe(:)*rc**(2*ic - 2)*cval
1450 404 : den(:, 2) = den(:, 2) - fe(:)*rc**(2*ic - 1)*coa
1451 6 : IF (ic > 1) THEN
1452 202 : den(:, 2) = den(:, 2) + REAL(2*ic - 2, dp)*fe(:)*rc**(2*ic - 3)*coa
1453 : END IF
1454 : END DO
1455 : END DO
1456 2 : DEALLOCATE (rc, fe)
1457 : END IF
1458 0 : ELSE IF (ASSOCIATED(sgp_potential)) THEN
1459 : CALL get_potential(potential=sgp_potential, has_nlcc=nlcc_type, &
1460 0 : n_nlcc=n_nlcc, a_nlcc=a_nlcc, c_nlcc=c_nlcc)
1461 0 : IF (nlcc_type) THEN
1462 0 : nr = qs_kind%grid_atom%nr
1463 0 : rr => qs_kind%grid_atom%rad
1464 0 : ALLOCATE (qs_kind%nlcc_pot(nr, 2), rc(nr), fe(nr))
1465 0 : den => qs_kind%nlcc_pot
1466 0 : den = 0.0_dp
1467 0 : DO i = 1, n_nlcc
1468 0 : alpha = a_nlcc(i)
1469 0 : fe(:) = EXP(-alpha*rr(:)*rr(:))
1470 0 : cval = c_nlcc(i)
1471 0 : den(:, 1) = den(:, 1) + cval*fe(:)
1472 0 : den(:, 2) = den(:, 2) - 2.0_dp*alpha*cval*rr(:)*fe(:)
1473 : END DO
1474 0 : DEALLOCATE (rc, fe)
1475 : END IF
1476 : ELSE
1477 : ! skip
1478 : END IF
1479 : END IF
1480 : END DO
1481 : END IF
1482 : ELSE
1483 0 : CPABORT("The pointer qs_kind_set is not associated")
1484 : END IF
1485 :
1486 1004 : END SUBROUTINE init_gapw_nlcc
1487 :
1488 : ! **************************************************************************************************
1489 : !> \brief Read an atomic kind data set from the input file.
1490 : !> \param qs_kind ...
1491 : !> \param kind_section ...
1492 : !> \param para_env ...
1493 : !> \param force_env_section ...
1494 : !> \param no_fail ...
1495 : !> \param method_id ...
1496 : !> \par History
1497 : !> - Creation (09.02.2002,MK)
1498 : !> - 20.09.2002,gt: adapted for POL/KG use (elp_potential)
1499 : !> - 05.03.2010: split elp_potential into fist_potential and kg_potential
1500 : ! **************************************************************************************************
1501 12823 : SUBROUTINE read_qs_kind(qs_kind, kind_section, para_env, force_env_section, no_fail, method_id)
1502 :
1503 : TYPE(qs_kind_type), INTENT(INOUT) :: qs_kind
1504 : TYPE(section_vals_type), POINTER :: kind_section
1505 : TYPE(mp_para_env_type), POINTER :: para_env
1506 : TYPE(section_vals_type), POINTER :: force_env_section
1507 : LOGICAL, INTENT(IN) :: no_fail
1508 : INTEGER, INTENT(IN) :: method_id
1509 :
1510 : CHARACTER(LEN=*), PARAMETER :: routineN = 'read_qs_kind'
1511 : INTEGER, PARAMETER :: maxbas = 20
1512 :
1513 : CHARACTER(LEN=2) :: element_symbol
1514 : CHARACTER(len=default_path_length) :: kg_potential_fn_kind, &
1515 : potential_file_name, potential_fn_kind
1516 : CHARACTER(LEN=default_string_length) :: akind_name, basis_type, keyword, &
1517 : kgpot_name, kgpot_type, &
1518 : potential_name, potential_type, tmp
1519 : CHARACTER(LEN=default_string_length), DIMENSION(4) :: description
1520 : CHARACTER(LEN=default_string_length), &
1521 12823 : DIMENSION(:), POINTER :: tmpstringlist
1522 : CHARACTER(LEN=default_string_length), &
1523 : DIMENSION(maxbas) :: basis_set_form, basis_set_name, &
1524 : basis_set_type
1525 : INTEGER :: handle, i, i_rep, iounit, ipaodesc, ipaopot, ipos, j, jj, k_rep, l, m, n_rep, &
1526 : nb_rep, nexp, ngauss, nlcc, nloc, nnl, norbitals, npaodesc, npaopot, nppnl, nspin, nu, z
1527 25646 : INTEGER, DIMENSION(:), POINTER :: add_el, elec_conf, orbitals
1528 : LOGICAL :: check, ecp_semi_local, explicit, explicit_basis, explicit_J, explicit_kgpot, &
1529 : explicit_potential, explicit_U, explicit_u_m_j, nobasis, section_enabled, &
1530 : subsection_enabled, update_input
1531 : REAL(KIND=dp) :: alpha, ccore, r, rc, zeff_correction
1532 : REAL(KIND=dp), DIMENSION(6) :: error
1533 25646 : REAL(KIND=dp), DIMENSION(:), POINTER :: a_nl, aloc, anlcc, cloc, cnlcc, nelec
1534 12823 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: h_nl
1535 12823 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: c_nl
1536 : TYPE(atom_ecppot_type) :: ecppot
1537 : TYPE(atom_sgp_potential_type) :: sgppot
1538 1346415 : TYPE(atom_upfpot_type) :: upfpot
1539 : TYPE(cp_logger_type), POINTER :: logger
1540 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set, sup_basis_set, &
1541 : tmp_basis_set
1542 : TYPE(section_vals_type), POINTER :: basis_section, bs_section, dft_plus_u_section, &
1543 : dft_section, enforce_occupation_section, kgpot_section, pao_desc_section, &
1544 : pao_pot_section, potential_section, spin_section
1545 : TYPE(sto_basis_set_type), POINTER :: sto_basis_set
1546 :
1547 12823 : CALL timeset(routineN, handle)
1548 :
1549 12823 : NULLIFY (logger)
1550 12823 : logger => cp_get_default_logger()
1551 12823 : iounit = cp_logger_get_default_io_unit(logger)
1552 :
1553 12823 : NULLIFY (elec_conf)
1554 :
1555 12823 : update_input = .TRUE.
1556 269283 : basis_set_name(:) = ""
1557 269283 : basis_set_type(:) = ""
1558 269283 : basis_set_form(:) = ""
1559 12823 : potential_name = ""
1560 12823 : potential_type = ""
1561 12823 : kgpot_name = ""
1562 12823 : kgpot_type = ""
1563 12823 : z = -1
1564 12823 : zeff_correction = 0.0_dp
1565 12823 : explicit = .FALSE.
1566 12823 : explicit_basis = .FALSE.
1567 12823 : explicit_J = .FALSE.
1568 12823 : explicit_kgpot = .FALSE.
1569 12823 : explicit_potential = .FALSE.
1570 12823 : explicit_U = .FALSE.
1571 12823 : explicit_u_m_j = .FALSE.
1572 :
1573 12823 : dft_section => section_vals_get_subs_vals(force_env_section, "DFT")
1574 12823 : CALL section_vals_get(kind_section, n_repetition=n_rep)
1575 12823 : k_rep = -1
1576 12823 : akind_name = qs_kind%name
1577 12823 : CALL uppercase(akind_name)
1578 : ! First we use the atom_name to find out the proper KIND section
1579 19246 : DO i_rep = 1, n_rep
1580 : CALL section_vals_val_get(kind_section, "_SECTION_PARAMETERS_", &
1581 15870 : c_val=keyword, i_rep_section=i_rep)
1582 15870 : CALL uppercase(keyword)
1583 19246 : IF (keyword == akind_name) THEN
1584 9447 : k_rep = i_rep
1585 9447 : EXIT
1586 : END IF
1587 : END DO
1588 : ! The search for the KIND section failed.. check for a QM/MM link atom
1589 12823 : IF (k_rep < 1) THEN
1590 3376 : ipos = INDEX(qs_kind%name, "_")
1591 3376 : IF (((ipos == 2) .OR. (ipos == 3)) .AND. (INDEX(qs_kind%name, "_ghost") == 0)) THEN
1592 : ! If the atm_name could not match any KIND section it maybe be a QM/MM link atom.
1593 : ! ghost atoms will be treated differently.
1594 64 : akind_name = qs_kind%name(1:ipos - 1)
1595 64 : CALL uppercase(akind_name)
1596 64 : DO i_rep = 1, n_rep
1597 : CALL section_vals_val_get(kind_section, "_SECTION_PARAMETERS_", &
1598 52 : c_val=keyword, i_rep_section=i_rep)
1599 52 : CALL uppercase(keyword)
1600 64 : IF (keyword == akind_name) THEN
1601 52 : k_rep = i_rep
1602 52 : EXIT
1603 : END IF
1604 : END DO
1605 : END IF
1606 : END IF
1607 : ! The search for the KIND section failed.. check element_symbol
1608 12823 : IF (k_rep < 1) THEN
1609 : ! If it's not a link atom let's check for the element and map
1610 : ! the KIND section to the element.
1611 3324 : element_symbol = qs_kind%element_symbol(1:2)
1612 3324 : CALL uppercase(element_symbol)
1613 3412 : DO i_rep = 1, n_rep
1614 : CALL section_vals_val_get(kind_section, "_SECTION_PARAMETERS_", &
1615 108 : c_val=keyword, i_rep_section=i_rep)
1616 108 : CALL uppercase(keyword)
1617 3412 : IF (keyword == element_symbol) THEN
1618 20 : k_rep = i_rep
1619 20 : EXIT
1620 : END IF
1621 : END DO
1622 : END IF
1623 : ! In case it should not really match any possible KIND section
1624 : ! let's look if a default one is defined..
1625 12823 : IF (k_rep < 1) THEN
1626 3320 : DO i_rep = 1, n_rep
1627 : CALL section_vals_val_get(kind_section, "_SECTION_PARAMETERS_", &
1628 68 : c_val=keyword, i_rep_section=i_rep)
1629 68 : CALL uppercase(keyword)
1630 3320 : IF (keyword == "DEFAULT") THEN
1631 52 : update_input = .FALSE.
1632 52 : k_rep = i_rep
1633 52 : EXIT
1634 : END IF
1635 : END DO
1636 : END IF
1637 12823 : IF (k_rep < 0 .AND. (.NOT. no_fail)) THEN
1638 : CALL cp_abort(__LOCATION__, &
1639 : "No &KIND section was possible to associate to the atomic kind <"// &
1640 : TRIM(akind_name)//">. The KIND section were also scanned for the"// &
1641 : " corresponding element <"//TRIM(qs_kind%element_symbol)//">"// &
1642 0 : " and for the DEFAULT section but no match was found. Check your input file!")
1643 : END IF
1644 : ! Retrieve information on element
1645 12823 : CALL get_ptable_info(qs_kind%element_symbol, ielement=z)
1646 :
1647 : ! Normal parsing of the KIND section
1648 12823 : IF (k_rep > 0) THEN
1649 : ! new style basis set input
1650 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1651 : keyword_name="BASIS_SET", &
1652 : explicit=explicit, &
1653 9571 : n_rep_val=nb_rep)
1654 9571 : IF (.NOT. explicit) nb_rep = 0
1655 9571 : CPASSERT(nb_rep <= maxbas)
1656 20939 : DO i = 1, nb_rep
1657 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1658 11368 : keyword_name="BASIS_SET", i_rep_val=i, c_vals=tmpstringlist)
1659 11368 : IF (SIZE(tmpstringlist) == 1) THEN
1660 : ! default is orbital type and GTO
1661 8613 : basis_set_type(i) = "ORB"
1662 8613 : basis_set_form(i) = "GTO"
1663 8613 : basis_set_name(i) = tmpstringlist(1)
1664 2755 : ELSEIF (SIZE(tmpstringlist) == 2) THEN
1665 : ! default is GTO
1666 2751 : basis_set_type(i) = tmpstringlist(1)
1667 2751 : basis_set_form(i) = "GTO"
1668 2751 : basis_set_name(i) = tmpstringlist(2)
1669 4 : ELSEIF (SIZE(tmpstringlist) == 3) THEN
1670 4 : basis_set_type(i) = tmpstringlist(1)
1671 4 : basis_set_form(i) = tmpstringlist(2)
1672 4 : basis_set_name(i) = tmpstringlist(3)
1673 : ELSE
1674 : CALL cp_abort(__LOCATION__, &
1675 0 : "invalid number of BASIS_SET keyword parameters: BASIS_SET [<TYPE>] [<FORM>] <NAME>")
1676 : END IF
1677 : ! check that we have a valid basis set form
1678 20939 : IF (basis_set_form(i) /= "GTO" .AND. basis_set_form(i) /= "STO") THEN
1679 0 : CPABORT("invalid BASIS_SET FORM parameter")
1680 : END IF
1681 : END DO
1682 :
1683 : ! parse PAO keywords
1684 : CALL section_vals_val_get(kind_section, keyword_name="PAO_BASIS_SIZE", i_rep_section=k_rep, &
1685 9571 : i_val=qs_kind%pao_basis_size)
1686 : CALL section_vals_val_get(kind_section, keyword_name="PAO_MODEL_FILE", i_rep_section=k_rep, &
1687 9571 : explicit=explicit)
1688 9571 : IF (explicit) THEN
1689 : CALL section_vals_val_get(kind_section, keyword_name="PAO_MODEL_FILE", i_rep_section=k_rep, &
1690 4 : c_val=qs_kind%pao_model_file)
1691 : END IF
1692 :
1693 : ! parse PAO_POTENTIAL sections
1694 9571 : pao_pot_section => section_vals_get_subs_vals(kind_section, "PAO_POTENTIAL", i_rep_section=k_rep)
1695 9571 : CALL section_vals_get(pao_pot_section, n_repetition=npaopot)
1696 19264 : ALLOCATE (qs_kind%pao_potentials(npaopot))
1697 9633 : DO ipaopot = 1, npaopot
1698 : CALL section_vals_val_get(pao_pot_section, keyword_name="MAXL", i_rep_section=ipaopot, &
1699 62 : i_val=qs_kind%pao_potentials(ipaopot)%maxl)
1700 : CALL section_vals_val_get(pao_pot_section, keyword_name="MAX_PROJECTOR", i_rep_section=ipaopot, &
1701 62 : i_val=qs_kind%pao_potentials(ipaopot)%max_projector)
1702 : CALL section_vals_val_get(pao_pot_section, keyword_name="BETA", i_rep_section=ipaopot, &
1703 62 : r_val=qs_kind%pao_potentials(ipaopot)%beta)
1704 : CALL section_vals_val_get(pao_pot_section, keyword_name="WEIGHT", i_rep_section=ipaopot, &
1705 9633 : r_val=qs_kind%pao_potentials(ipaopot)%weight)
1706 : END DO
1707 :
1708 : ! parse PAO_DESCRIPTOR sections
1709 9571 : pao_desc_section => section_vals_get_subs_vals(kind_section, "PAO_DESCRIPTOR", i_rep_section=k_rep)
1710 9571 : CALL section_vals_get(pao_desc_section, n_repetition=npaodesc)
1711 19172 : ALLOCATE (qs_kind%pao_descriptors(npaodesc))
1712 9589 : DO ipaodesc = 1, npaodesc
1713 : CALL section_vals_val_get(pao_desc_section, keyword_name="BETA", i_rep_section=ipaodesc, &
1714 18 : r_val=qs_kind%pao_descriptors(ipaodesc)%beta)
1715 : CALL section_vals_val_get(pao_desc_section, keyword_name="SCREENING", i_rep_section=ipaodesc, &
1716 18 : r_val=qs_kind%pao_descriptors(ipaodesc)%screening)
1717 : CALL section_vals_val_get(pao_desc_section, keyword_name="WEIGHT", i_rep_section=ipaodesc, &
1718 9589 : r_val=qs_kind%pao_descriptors(ipaodesc)%weight)
1719 : END DO
1720 :
1721 : ! parse ELEC_CONF
1722 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1723 9571 : keyword_name="ELEC_CONF", n_rep_val=i)
1724 9571 : IF (i > 0) THEN
1725 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1726 4 : keyword_name="ELEC_CONF", i_vals=elec_conf)
1727 4 : CALL set_qs_kind(qs_kind, elec_conf=elec_conf)
1728 : END IF
1729 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1730 9571 : keyword_name="CORE_CORRECTION", r_val=zeff_correction)
1731 : ! parse POTENTIAL
1732 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1733 9571 : keyword_name="POTENTIAL_FILE_NAME", c_val=potential_fn_kind)
1734 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1735 9571 : keyword_name="POTENTIAL_TYPE", c_val=potential_type)
1736 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1737 9571 : explicit=explicit, keyword_name="POTENTIAL", c_vals=tmpstringlist)
1738 9571 : IF (explicit) THEN
1739 9303 : IF (SIZE(tmpstringlist) == 1) THEN
1740 : ! old type of input: start of name defines type
1741 9259 : potential_name = tmpstringlist(1)
1742 9259 : IF (potential_type == "") THEN
1743 9259 : ipos = INDEX(potential_name, "-")
1744 9259 : IF (ipos > 1) THEN
1745 8231 : potential_type = potential_name(:ipos - 1)
1746 : ELSE
1747 1028 : potential_type = potential_name
1748 : END IF
1749 : END IF
1750 44 : ELSEIF (SIZE(tmpstringlist) == 2) THEN
1751 44 : potential_type = tmpstringlist(1)
1752 44 : potential_name = tmpstringlist(2)
1753 : ELSE
1754 0 : CPABORT("POTENTIAL input list is not correct")
1755 : END IF
1756 : END IF
1757 9571 : CALL uppercase(potential_type)
1758 :
1759 : ! Parse KG POTENTIAL
1760 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1761 9571 : keyword_name="KG_POTENTIAL_FILE_NAME", c_val=kg_potential_fn_kind)
1762 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1763 9571 : keyword_name="KG_POTENTIAL", c_val=kgpot_name)
1764 :
1765 : ! Semi-local vs. full nonlocal form of ECPs
1766 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1767 9571 : keyword_name="ECP_SEMI_LOCAL", l_val=ecp_semi_local)
1768 :
1769 : ! Assign atomic covalent radius
1770 9571 : qs_kind%covalent_radius = ptable(z)%covalent_radius*bohr
1771 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1772 9571 : keyword_name="COVALENT_RADIUS", r_val=r)
1773 9571 : IF (r > 0.0_dp) qs_kind%covalent_radius = r
1774 :
1775 : ! Assign atomic van der Waals radius
1776 9571 : qs_kind%vdw_radius = ptable(z)%vdw_radius*bohr
1777 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1778 9571 : keyword_name="VDW_RADIUS", r_val=r)
1779 9571 : IF (r > 0.0_dp) qs_kind%vdw_radius = r
1780 :
1781 : ! Assign atom dependent defaults, only H special case
1782 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, n_rep_val=i, &
1783 9571 : keyword_name="HARD_EXP_RADIUS")
1784 9571 : IF (i == 0) THEN
1785 9517 : IF (z == 1) THEN
1786 4168 : qs_kind%hard_radius = 1.2_dp
1787 : ELSE
1788 5349 : qs_kind%hard_radius = 0.8_dp*bohr
1789 : END IF
1790 : ELSE
1791 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1792 54 : keyword_name="HARD_EXP_RADIUS", r_val=qs_kind%hard_radius)
1793 : END IF
1794 :
1795 : ! assign atom dependent defaults, only H special case
1796 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, n_rep_val=i, &
1797 9571 : keyword_name="RHO0_EXP_RADIUS")
1798 9571 : IF (i == 0) THEN
1799 9571 : qs_kind%hard0_radius = qs_kind%hard_radius
1800 : ELSE
1801 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1802 0 : keyword_name="RHO0_EXP_RADIUS", r_val=qs_kind%hard0_radius)
1803 : END IF
1804 9571 : IF (qs_kind%hard_radius < qs_kind%hard0_radius) &
1805 0 : CPABORT("rc0 should be <= rc")
1806 :
1807 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1808 9571 : keyword_name="MAX_RAD_LOCAL", r_val=qs_kind%max_rad_local)
1809 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1810 9571 : keyword_name="LEBEDEV_GRID", i_val=qs_kind%ngrid_ang)
1811 9571 : IF (qs_kind%ngrid_ang <= 0) &
1812 0 : CPABORT("# point lebedev grid < 0")
1813 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1814 9571 : keyword_name="RADIAL_GRID", i_val=qs_kind%ngrid_rad)
1815 9571 : IF (qs_kind%ngrid_rad <= 0) &
1816 0 : CPABORT("# point radial grid < 0")
1817 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1818 9571 : keyword_name="GPW_TYPE", l_val=qs_kind%gpw_type_forced)
1819 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1820 9571 : keyword_name="GHOST", l_val=qs_kind%ghost)
1821 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1822 9571 : keyword_name="FLOATING_BASIS_CENTER", l_val=qs_kind%floating)
1823 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1824 9571 : keyword_name="NO_OPTIMIZE", l_val=qs_kind%no_optimize)
1825 :
1826 : ! Magnetization
1827 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1828 9571 : keyword_name="MAGNETIZATION", r_val=qs_kind%magnetization)
1829 : ! DFTB3 param
1830 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1831 9571 : keyword_name="DFTB3_PARAM", r_val=qs_kind%dudq_dftb3)
1832 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1833 9571 : keyword_name="LMAX_DFTB", i_val=qs_kind%lmax_dftb)
1834 :
1835 : ! MAOS
1836 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1837 9571 : keyword_name="MAO", i_val=qs_kind%mao)
1838 :
1839 : ! Read the BS subsection of the current atomic kind, if enabled
1840 9571 : NULLIFY (bs_section)
1841 : bs_section => section_vals_get_subs_vals(kind_section, "BS", &
1842 9571 : i_rep_section=k_rep)
1843 : section_enabled = .FALSE.
1844 : CALL section_vals_val_get(bs_section, "_SECTION_PARAMETERS_", &
1845 9571 : l_val=section_enabled)
1846 9571 : IF (section_enabled) THEN
1847 : ! test for conflict with magnetization
1848 60 : IF (qs_kind%magnetization /= 0.0_dp) THEN
1849 : CALL cp_abort(__LOCATION__, "BS Section is in conflict with non-zero magnetization "// &
1850 0 : "for this atom kind.")
1851 : END IF
1852 60 : qs_kind%bs_occupation = .TRUE.
1853 : !Alpha spin
1854 60 : NULLIFY (spin_section)
1855 60 : spin_section => section_vals_get_subs_vals(bs_section, "ALPHA")
1856 60 : CALL section_vals_get(spin_section, explicit=explicit)
1857 60 : IF (explicit) THEN
1858 60 : NULLIFY (add_el)
1859 : CALL section_vals_val_get(spin_section, &
1860 60 : keyword_name="NEL", i_vals=add_el)
1861 60 : CPASSERT(ASSOCIATED(add_el))
1862 180 : ALLOCATE (qs_kind%addel(SIZE(add_el), 2))
1863 328 : qs_kind%addel = 0
1864 134 : qs_kind%addel(1:SIZE(add_el), 1) = add_el(1:SIZE(add_el))
1865 60 : NULLIFY (add_el)
1866 : CALL section_vals_val_get(spin_section, &
1867 60 : keyword_name="L", i_vals=add_el)
1868 60 : CPASSERT(ASSOCIATED(add_el))
1869 60 : CPASSERT(SIZE(add_el) == SIZE(qs_kind%addel, 1))
1870 180 : ALLOCATE (qs_kind%laddel(SIZE(add_el), 2))
1871 328 : qs_kind%laddel = 0
1872 134 : qs_kind%laddel(1:SIZE(add_el), 1) = add_el(1:SIZE(add_el))
1873 180 : ALLOCATE (qs_kind%naddel(SIZE(add_el), 2))
1874 328 : qs_kind%naddel = 0
1875 60 : NULLIFY (add_el)
1876 : CALL section_vals_val_get(spin_section, &
1877 60 : keyword_name="N", n_rep_val=i)
1878 60 : IF (i > 0) THEN
1879 : CALL section_vals_val_get(spin_section, &
1880 60 : keyword_name="N", i_vals=add_el)
1881 60 : IF (SIZE(add_el) == SIZE(qs_kind%addel, 1)) THEN
1882 134 : qs_kind%naddel(1:SIZE(add_el), 1) = add_el(1:SIZE(add_el))
1883 : END IF
1884 : END IF
1885 : END IF
1886 : ! Beta spin
1887 60 : NULLIFY (spin_section)
1888 60 : spin_section => section_vals_get_subs_vals(bs_section, "BETA")
1889 60 : CALL section_vals_get(spin_section, explicit=explicit)
1890 60 : IF (explicit) THEN
1891 60 : NULLIFY (add_el)
1892 : CALL section_vals_val_get(spin_section, &
1893 60 : keyword_name="NEL", i_vals=add_el)
1894 60 : CPASSERT(SIZE(add_el) == SIZE(qs_kind%addel, 1))
1895 134 : qs_kind%addel(1:SIZE(add_el), 2) = add_el(1:SIZE(add_el))
1896 328 : qs_kind%addel(:, :) = qs_kind%addel(:, :)
1897 60 : NULLIFY (add_el)
1898 : CALL section_vals_val_get(spin_section, &
1899 60 : keyword_name="L", i_vals=add_el)
1900 60 : CPASSERT(SIZE(add_el) == SIZE(qs_kind%addel, 1))
1901 134 : qs_kind%laddel(1:SIZE(add_el), 2) = add_el(1:SIZE(add_el))
1902 :
1903 : CALL section_vals_val_get(spin_section, &
1904 60 : keyword_name="N", n_rep_val=i)
1905 60 : IF (i > 0) THEN
1906 60 : NULLIFY (add_el)
1907 : CALL section_vals_val_get(spin_section, &
1908 60 : keyword_name="N", i_vals=add_el)
1909 60 : IF (SIZE(add_el) == SIZE(qs_kind%addel, 1)) THEN
1910 134 : qs_kind%naddel(1:SIZE(add_el), 2) = add_el(1:SIZE(add_el))
1911 : END IF
1912 : END IF
1913 : END IF
1914 : END IF
1915 :
1916 : ! Read the DFT+U subsection of the current atomic kind, if enabled
1917 :
1918 9571 : NULLIFY (dft_plus_u_section)
1919 : dft_plus_u_section => section_vals_get_subs_vals(kind_section, &
1920 : subsection_name="DFT_PLUS_U", &
1921 9571 : i_rep_section=k_rep)
1922 : section_enabled = .FALSE.
1923 : CALL section_vals_val_get(dft_plus_u_section, &
1924 : keyword_name="_SECTION_PARAMETERS_", &
1925 9571 : l_val=section_enabled)
1926 124423 : IF (section_enabled) THEN
1927 32 : ALLOCATE (qs_kind%dft_plus_u)
1928 : NULLIFY (qs_kind%dft_plus_u%nelec)
1929 : NULLIFY (qs_kind%dft_plus_u%orbitals)
1930 : CALL section_vals_val_get(dft_plus_u_section, &
1931 : keyword_name="L", &
1932 32 : i_val=l)
1933 32 : qs_kind%dft_plus_u%l = l
1934 : #if defined(__SIRIUS)
1935 : CALL section_vals_val_get(dft_plus_u_section, &
1936 : keyword_name="N", &
1937 32 : i_val=nu)
1938 32 : qs_kind%dft_plus_u%n = nu
1939 :
1940 : CALL section_vals_val_get(dft_plus_u_section, &
1941 : keyword_name="U", &
1942 : r_val=qs_kind%dft_plus_u%U, &
1943 32 : explicit=explicit_U)
1944 :
1945 : CALL section_vals_val_get(dft_plus_u_section, &
1946 : keyword_name="J", &
1947 : r_val=qs_kind%dft_plus_u%J, &
1948 32 : explicit=explicit_J)
1949 :
1950 : CALL section_vals_val_get(dft_plus_u_section, &
1951 : keyword_name="alpha", &
1952 32 : r_val=qs_kind%dft_plus_u%alpha)
1953 :
1954 : CALL section_vals_val_get(dft_plus_u_section, &
1955 : keyword_name="beta", &
1956 32 : r_val=qs_kind%dft_plus_u%beta)
1957 :
1958 : CALL section_vals_val_get(dft_plus_u_section, &
1959 : keyword_name="J0", &
1960 32 : r_val=qs_kind%dft_plus_u%J0)
1961 :
1962 : CALL section_vals_val_get(dft_plus_u_section, &
1963 : keyword_name="occupation", &
1964 32 : r_val=qs_kind%dft_plus_u%occupation)
1965 : #else
1966 : nu = 0
1967 : #endif
1968 :
1969 : CALL section_vals_val_get(dft_plus_u_section, &
1970 : keyword_name="U_MINUS_J", &
1971 : r_val=qs_kind%dft_plus_u%u_minus_j_target, &
1972 32 : explicit=explicit_u_m_j)
1973 :
1974 32 : IF ((explicit_U .OR. explicit_J) .AND. explicit_u_m_j) THEN
1975 0 : CPABORT("DFT+U| specifying U or J and U_MINUS_J parameters are mutually exclusive.")
1976 : END IF
1977 :
1978 : CALL section_vals_val_get(dft_plus_u_section, &
1979 : keyword_name="U_RAMPING", &
1980 32 : r_val=qs_kind%dft_plus_u%u_ramping)
1981 : CALL section_vals_val_get(dft_plus_u_section, &
1982 : keyword_name="INIT_U_RAMPING_EACH_SCF", &
1983 32 : l_val=qs_kind%dft_plus_u%init_u_ramping_each_scf)
1984 32 : IF (qs_kind%dft_plus_u%u_ramping > 0.0_dp) THEN
1985 8 : qs_kind%dft_plus_u%u_minus_j = 0.0_dp
1986 : ELSE
1987 24 : qs_kind%dft_plus_u%u_minus_j = qs_kind%dft_plus_u%u_minus_j_target
1988 : END IF
1989 : CALL section_vals_val_get(dft_plus_u_section, &
1990 : keyword_name="EPS_U_RAMPING", &
1991 32 : r_val=qs_kind%dft_plus_u%eps_u_ramping)
1992 :
1993 32 : NULLIFY (enforce_occupation_section)
1994 : enforce_occupation_section => section_vals_get_subs_vals(dft_plus_u_section, &
1995 32 : subsection_name="ENFORCE_OCCUPATION")
1996 : subsection_enabled = .FALSE.
1997 : CALL section_vals_val_get(enforce_occupation_section, &
1998 : keyword_name="_SECTION_PARAMETERS_", &
1999 32 : l_val=subsection_enabled)
2000 32 : IF (subsection_enabled) THEN
2001 4 : NULLIFY (nelec)
2002 : CALL section_vals_val_get(enforce_occupation_section, &
2003 : keyword_name="NELEC", &
2004 4 : r_vals=nelec)
2005 4 : nspin = SIZE(nelec)
2006 12 : ALLOCATE (qs_kind%dft_plus_u%nelec(nspin))
2007 8 : qs_kind%dft_plus_u%nelec(:) = nelec(:)
2008 4 : NULLIFY (orbitals)
2009 : CALL section_vals_val_get(enforce_occupation_section, &
2010 : keyword_name="ORBITALS", &
2011 4 : i_vals=orbitals)
2012 4 : norbitals = SIZE(orbitals)
2013 4 : IF (norbitals <= 0 .OR. norbitals > 2*l + 1) &
2014 : CALL cp_abort(__LOCATION__, "DFT+U| Invalid number of ORBITALS specified: "// &
2015 0 : "1 to 2*L+1 integer numbers are expected")
2016 12 : ALLOCATE (qs_kind%dft_plus_u%orbitals(norbitals))
2017 16 : qs_kind%dft_plus_u%orbitals(:) = orbitals(:)
2018 4 : NULLIFY (orbitals)
2019 16 : DO m = 1, norbitals
2020 12 : IF (qs_kind%dft_plus_u%orbitals(m) > l) &
2021 0 : CPABORT("DFT+U| Invalid orbital magnetic quantum number specified: m > l")
2022 12 : IF (qs_kind%dft_plus_u%orbitals(m) < -l) &
2023 0 : CPABORT("DFT+U| Invalid orbital magnetic quantum number specified: m < -l")
2024 52 : DO j = 1, norbitals
2025 48 : IF (j /= m) THEN
2026 24 : IF (qs_kind%dft_plus_u%orbitals(j) == qs_kind%dft_plus_u%orbitals(m)) &
2027 0 : CPABORT("DFT+U| An orbital magnetic quantum number was specified twice")
2028 : END IF
2029 : END DO
2030 : END DO
2031 : CALL section_vals_val_get(enforce_occupation_section, &
2032 : keyword_name="EPS_SCF", &
2033 4 : r_val=qs_kind%dft_plus_u%eps_scf)
2034 : CALL section_vals_val_get(enforce_occupation_section, &
2035 : keyword_name="MAX_SCF", &
2036 4 : i_val=i)
2037 4 : qs_kind%dft_plus_u%max_scf = MAX(-1, i)
2038 : CALL section_vals_val_get(enforce_occupation_section, &
2039 : keyword_name="SMEAR", &
2040 4 : l_val=qs_kind%dft_plus_u%smear)
2041 : END IF ! subsection enabled
2042 : END IF ! section enabled
2043 :
2044 : END IF
2045 :
2046 : ! Allocate and initialise the orbital basis set data set structure
2047 12823 : CALL init_orbital_pointers(5) ! debug the SUN optimizer
2048 :
2049 : ! BASIS and POTENTIAL read only when strictly necessary otherwise, even if not used
2050 : ! we just print misleading informations
2051 12823 : explicit_basis = .FALSE.
2052 12823 : IF (k_rep > 0) THEN
2053 : basis_section => section_vals_get_subs_vals(kind_section, "BASIS", i_rep_section=k_rep, &
2054 9571 : can_return_null=.TRUE.)
2055 9571 : CALL section_vals_get(basis_section, explicit=explicit_basis)
2056 : END IF
2057 :
2058 12823 : explicit_potential = .FALSE.
2059 12823 : IF (k_rep > 0) THEN
2060 : potential_section => section_vals_get_subs_vals(kind_section, "POTENTIAL", &
2061 9571 : i_rep_section=k_rep, can_return_null=.TRUE.)
2062 9571 : CALL section_vals_get(potential_section, explicit=explicit_potential)
2063 : END IF
2064 :
2065 12823 : explicit_kgpot = .FALSE.
2066 12823 : IF (k_rep > 0) THEN
2067 : kgpot_section => section_vals_get_subs_vals(kind_section, "KG_POTENTIAL", &
2068 9571 : i_rep_section=k_rep, can_return_null=.TRUE.)
2069 9571 : CALL section_vals_get(kgpot_section, explicit=explicit_kgpot)
2070 : END IF
2071 :
2072 15063 : SELECT CASE (method_id)
2073 : CASE (do_method_rm1, do_method_am1, do_method_mndo, do_method_pdg, do_method_pm3, do_method_pm6, &
2074 : do_method_pm6fm, do_method_mndod, do_method_pnnl)
2075 : ! Allocate all_potential
2076 2240 : CALL allocate_potential(qs_kind%all_potential)
2077 2240 : CALL set_default_all_potential(qs_kind%all_potential, z, zeff_correction)
2078 2240 : CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
2079 2240 : IF (.NOT. ASSOCIATED(elec_conf)) THEN
2080 2240 : CALL get_potential(potential=qs_kind%all_potential, elec_conf=elec_conf)
2081 2240 : CALL set_qs_kind(qs_kind, elec_conf=elec_conf)
2082 : END IF
2083 2240 : CPASSERT(.NOT. qs_kind%floating)
2084 2240 : IF (qs_kind%ghost) THEN
2085 0 : CALL get_qs_kind(qs_kind=qs_kind, elec_conf=elec_conf)
2086 0 : elec_conf(:) = 0
2087 : CALL get_potential(potential=qs_kind%all_potential, &
2088 0 : elec_conf=elec_conf)
2089 0 : elec_conf(:) = 0
2090 : CALL set_potential(potential=qs_kind%all_potential, &
2091 : zeff=0.0_dp, &
2092 0 : zeff_correction=0.0_dp)
2093 : END IF
2094 :
2095 : ! Basis set (Parameters)
2096 : ! Setup proper semiempirical parameters
2097 2240 : check = .NOT. ASSOCIATED(qs_kind%se_parameter)
2098 2240 : CPASSERT(check)
2099 2240 : CALL semi_empirical_create(qs_kind%se_parameter)
2100 : ! Check if we allow p-orbitals on H
2101 438 : SELECT CASE (z)
2102 : CASE (1)
2103 2240 : IF (k_rep > 0) THEN
2104 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
2105 52 : keyword_name="SE_P_ORBITALS_ON_H", l_val=qs_kind%se_parameter%p_orbitals_on_h)
2106 : END IF
2107 : CASE DEFAULT
2108 : ! No special cases for other elements..
2109 : END SELECT
2110 : ! Set default parameters
2111 2240 : CALL section_vals_val_get(dft_section, "QS%SE%STO_NG", i_val=ngauss)
2112 2240 : CALL se_param_set_default(qs_kind%se_parameter, z, method_id)
2113 2240 : NULLIFY (tmp_basis_set)
2114 2240 : CALL init_se_param(qs_kind%se_parameter, tmp_basis_set, ngauss)
2115 2240 : CALL add_basis_set_to_container(qs_kind%basis_sets, tmp_basis_set, "ORB")
2116 : CALL init_potential(qs_kind%all_potential, itype="BARE", &
2117 2240 : zeff=qs_kind%se_parameter%zeff, zeff_correction=zeff_correction)
2118 2240 : qs_kind%se_parameter%zeff = qs_kind%se_parameter%zeff - zeff_correction
2119 :
2120 2240 : check = ((potential_name /= '') .OR. explicit_potential)
2121 : IF (check) &
2122 : CALL cp_warn(__LOCATION__, &
2123 : "Information provided in the input file regarding POTENTIAL for KIND <"// &
2124 80 : TRIM(qs_kind%name)//"> will be ignored!")
2125 :
2126 2240 : check = ((k_rep > 0) .OR. explicit_basis)
2127 : IF (check) &
2128 : CALL cp_warn(__LOCATION__, &
2129 : "Information provided in the input file regarding BASIS for KIND <"// &
2130 116 : TRIM(qs_kind%name)//"> will be ignored!")
2131 :
2132 : CASE (do_method_dftb)
2133 : ! Allocate all_potential
2134 480 : CALL allocate_potential(qs_kind%all_potential)
2135 480 : CALL set_default_all_potential(qs_kind%all_potential, z, zeff_correction)
2136 480 : CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
2137 480 : IF (.NOT. ASSOCIATED(elec_conf)) THEN
2138 480 : CALL get_potential(potential=qs_kind%all_potential, elec_conf=elec_conf)
2139 480 : CALL set_qs_kind(qs_kind, elec_conf=elec_conf)
2140 : END IF
2141 480 : CPASSERT(.NOT. qs_kind%floating)
2142 480 : IF (qs_kind%ghost) THEN
2143 0 : CALL get_qs_kind(qs_kind=qs_kind, elec_conf=elec_conf)
2144 0 : elec_conf(:) = 0
2145 : CALL get_potential(potential=qs_kind%all_potential, &
2146 0 : elec_conf=elec_conf)
2147 0 : elec_conf(:) = 0
2148 : CALL set_potential(potential=qs_kind%all_potential, &
2149 : zeff=0.0_dp, &
2150 0 : zeff_correction=0.0_dp)
2151 : END IF
2152 :
2153 480 : check = ((potential_name /= '') .OR. explicit_potential)
2154 : IF (check) &
2155 : CALL cp_warn(__LOCATION__, &
2156 : "Information provided in the input file regarding POTENTIAL for KIND <"// &
2157 0 : TRIM(qs_kind%name)//"> will be ignored!")
2158 :
2159 480 : check = ((k_rep > 0) .OR. explicit_basis)
2160 : IF (check) &
2161 : CALL cp_warn(__LOCATION__, &
2162 : "Information provided in the input file regarding BASIS for KIND <"// &
2163 44 : TRIM(qs_kind%name)//"> will be ignored!")
2164 :
2165 : CASE (do_method_xtb)
2166 : ! Allocate all_potential
2167 692 : CALL allocate_potential(qs_kind%all_potential)
2168 692 : CALL set_default_all_potential(qs_kind%all_potential, z, zeff_correction)
2169 692 : CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
2170 692 : IF (.NOT. ASSOCIATED(elec_conf)) THEN
2171 692 : CALL get_potential(potential=qs_kind%all_potential, elec_conf=elec_conf)
2172 692 : CALL set_qs_kind(qs_kind, elec_conf=elec_conf)
2173 : END IF
2174 692 : CPASSERT(.NOT. qs_kind%floating)
2175 692 : IF (qs_kind%ghost) THEN
2176 0 : CALL get_qs_kind(qs_kind=qs_kind, elec_conf=elec_conf)
2177 0 : elec_conf(:) = 0
2178 : CALL get_potential(potential=qs_kind%all_potential, &
2179 0 : elec_conf=elec_conf)
2180 0 : elec_conf(:) = 0
2181 : CALL set_potential(potential=qs_kind%all_potential, &
2182 : zeff=0.0_dp, &
2183 0 : zeff_correction=0.0_dp)
2184 : END IF
2185 :
2186 692 : check = ((potential_name /= '') .OR. explicit_potential)
2187 : IF (check) &
2188 : CALL cp_warn(__LOCATION__, &
2189 : "Information provided in the input file regarding POTENTIAL for KIND <"// &
2190 0 : TRIM(qs_kind%name)//"> will be ignored!")
2191 :
2192 692 : check = ((k_rep > 0) .OR. explicit_basis)
2193 : IF (check) &
2194 : CALL cp_warn(__LOCATION__, &
2195 : "Information provided in the input file regarding BASIS for KIND <"// &
2196 0 : TRIM(qs_kind%name)//"> will be ignored!")
2197 :
2198 : CASE (do_method_pw)
2199 : ! PW DFT
2200 : ! Allocate and initialise the potential data set structure
2201 22 : IF (potential_name /= '') THEN
2202 22 : SELECT CASE (TRIM(potential_type))
2203 : CASE ("ALL", "ECP")
2204 : CALL cp_abort(__LOCATION__, &
2205 : "PW DFT calculations only with potential type UPF or GTH possible."// &
2206 : " <"//TRIM(potential_type)//"> was specified "// &
2207 0 : "for the atomic kind <"//TRIM(qs_kind%name))
2208 : CASE ("GTH")
2209 2 : IF (potential_fn_kind == "-") THEN
2210 2 : CALL section_vals_val_get(dft_section, "POTENTIAL_FILE_NAME", c_val=potential_file_name)
2211 : ELSE
2212 0 : potential_file_name = potential_fn_kind
2213 : END IF
2214 2 : CALL allocate_potential(qs_kind%gth_potential)
2215 : CALL read_potential(qs_kind%element_symbol, potential_name, &
2216 : qs_kind%gth_potential, zeff_correction, para_env, &
2217 2 : potential_file_name, potential_section, update_input)
2218 2 : CALL set_potential(qs_kind%gth_potential, z=z)
2219 2 : CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
2220 2 : IF (.NOT. ASSOCIATED(elec_conf)) THEN
2221 2 : CALL get_potential(potential=qs_kind%gth_potential, elec_conf=elec_conf)
2222 2 : CALL set_qs_kind(qs_kind, elec_conf=elec_conf)
2223 : ELSE
2224 0 : CALL set_potential(potential=qs_kind%gth_potential, elec_conf=elec_conf)
2225 : END IF
2226 : CASE ("UPF")
2227 2100 : ALLOCATE (qs_kind%upf_potential)
2228 20 : qs_kind%upf_potential%zion = 0
2229 20 : qs_kind%upf_potential%filename = ADJUSTL(TRIM(potential_name))
2230 20 : CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
2231 20 : IF (.NOT. ASSOCIATED(elec_conf)) THEN
2232 20 : CALL set_qs_kind(qs_kind, elec_conf=qs_kind%upf_potential%econf)
2233 : END IF
2234 : CASE DEFAULT
2235 : CALL cp_abort(__LOCATION__, &
2236 : "An invalid potential type <"// &
2237 : TRIM(potential_type)//"> was specified "// &
2238 : "for the atomic kind <"// &
2239 22 : TRIM(qs_kind%name))
2240 : END SELECT
2241 : ELSE
2242 : CALL cp_abort(__LOCATION__, &
2243 : "No potential type was defined for the "// &
2244 0 : "atomic kind <"//TRIM(qs_kind%name)//">")
2245 : END IF
2246 :
2247 : CASE DEFAULT
2248 :
2249 : ! set ngauss for STO expansion
2250 9389 : CALL section_vals_val_get(dft_section, "QS%STO_NG", i_val=ngauss)
2251 : ! Allocate and initialise the basis set data set structure
2252 : ! first external basis sets
2253 20675 : DO i = 1, nb_rep
2254 22568 : SELECT CASE (basis_set_form(i))
2255 : CASE ("GTO")
2256 11282 : NULLIFY (tmp_basis_set)
2257 11282 : CALL allocate_gto_basis_set(tmp_basis_set)
2258 : CALL read_gto_basis_set(qs_kind%element_symbol, basis_set_name(i), &
2259 11282 : tmp_basis_set, para_env, dft_section)
2260 : CASE ("STO")
2261 4 : NULLIFY (sto_basis_set)
2262 4 : CALL allocate_sto_basis_set(sto_basis_set)
2263 : CALL read_sto_basis_set(qs_kind%element_symbol, basis_set_name(i), &
2264 4 : sto_basis_set, para_env, dft_section)
2265 4 : NULLIFY (tmp_basis_set)
2266 4 : CALL create_gto_from_sto_basis(sto_basis_set, tmp_basis_set, ngauss)
2267 4 : CALL deallocate_sto_basis_set(sto_basis_set)
2268 : CASE DEFAULT
2269 : CALL cp_abort(__LOCATION__, &
2270 : "Invalid basis set form "//TRIM(basis_set_form(i))// &
2271 11286 : "for atomic kind <"//TRIM(qs_kind%name)//">")
2272 : END SELECT
2273 11286 : tmp = basis_set_type(i)
2274 11286 : CALL uppercase(tmp)
2275 20675 : CALL add_basis_set_to_container(qs_kind%basis_sets, tmp_basis_set, tmp)
2276 : END DO
2277 : ! now explicit basis sets
2278 9389 : IF (explicit_basis) THEN
2279 162 : CALL section_vals_get(basis_section, n_repetition=nexp)
2280 324 : DO i = 1, nexp
2281 162 : NULLIFY (tmp_basis_set)
2282 162 : CALL allocate_gto_basis_set(tmp_basis_set)
2283 : CALL read_gto_basis_set(qs_kind%element_symbol, basis_type, &
2284 162 : tmp_basis_set, basis_section, i, dft_section)
2285 162 : tmp = basis_type
2286 162 : CALL uppercase(tmp)
2287 324 : CALL add_basis_set_to_container(qs_kind%basis_sets, tmp_basis_set, tmp)
2288 : END DO
2289 : END IF
2290 : ! combine multiple basis sets
2291 197169 : DO i = 1, SIZE(qs_kind%basis_sets)
2292 187780 : NULLIFY (tmp_basis_set)
2293 : CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis_set, &
2294 187780 : inumbas=i, basis_type=basis_type)
2295 187780 : IF (basis_type == "") CYCLE
2296 11448 : jj = i
2297 226775 : DO j = i + 1, SIZE(qs_kind%basis_sets)
2298 215327 : jj = jj + 1
2299 215327 : NULLIFY (sup_basis_set)
2300 : CALL get_basis_from_container(qs_kind%basis_sets, basis_set=sup_basis_set, &
2301 215327 : inumbas=jj, basis_type=tmp)
2302 226775 : IF (basis_type == tmp) THEN
2303 : ! we found a match, combine the basis sets and delete the second
2304 0 : CALL combine_basis_sets(tmp_basis_set, sup_basis_set)
2305 0 : CALL remove_basis_from_container(qs_kind%basis_sets, jj)
2306 0 : jj = jj - 1
2307 : END IF
2308 : END DO
2309 197169 : NULLIFY (sup_basis_set)
2310 : END DO
2311 :
2312 : ! check that we have an orbital basis set
2313 9389 : nobasis = .TRUE.
2314 197169 : DO i = 1, SIZE(qs_kind%basis_sets)
2315 187780 : NULLIFY (tmp_basis_set)
2316 : CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis_set, &
2317 187780 : inumbas=i, basis_type=basis_type)
2318 197169 : IF (basis_type == "ORB") nobasis = .FALSE.
2319 : END DO
2320 9389 : IF (nobasis) THEN
2321 : CALL cp_abort(__LOCATION__, &
2322 : "No basis set type was defined for the "// &
2323 0 : "atomic kind <"//TRIM(qs_kind%name)//">")
2324 : END IF
2325 :
2326 : ! If Ghost atom we don't need to allocate/initialize anything connected to POTENTIAL
2327 24452 : IF (qs_kind%ghost .OR. qs_kind%floating) THEN
2328 142 : IF (ASSOCIATED(qs_kind%elec_conf)) qs_kind%elec_conf = 0
2329 : ELSE
2330 : ! Allocate and initialise the potential data set structure
2331 9247 : IF ((potential_name /= '') .OR. explicit_potential) THEN
2332 : ! determine the pseudopotential file to search
2333 9247 : IF (potential_fn_kind == "-") THEN
2334 9237 : CALL section_vals_val_get(dft_section, "POTENTIAL_FILE_NAME", c_val=potential_file_name)
2335 : ELSE
2336 10 : potential_file_name = potential_fn_kind
2337 : END IF
2338 : !
2339 10271 : SELECT CASE (TRIM(potential_type))
2340 : CASE ("ALL")
2341 1024 : CALL allocate_potential(qs_kind%all_potential)
2342 : CALL read_potential(qs_kind%element_symbol, potential_name, &
2343 : qs_kind%all_potential, zeff_correction, para_env, &
2344 1024 : potential_file_name, potential_section, update_input)
2345 1024 : CALL set_potential(qs_kind%all_potential, z=z)
2346 1024 : CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
2347 1024 : IF (.NOT. ASSOCIATED(elec_conf)) THEN
2348 1024 : CALL get_potential(potential=qs_kind%all_potential, elec_conf=elec_conf)
2349 1024 : CALL set_qs_kind(qs_kind, elec_conf=elec_conf)
2350 : ELSE
2351 0 : CALL set_potential(potential=qs_kind%all_potential, elec_conf=elec_conf)
2352 : END IF
2353 : CASE ("GTH")
2354 8199 : CALL allocate_potential(qs_kind%gth_potential)
2355 : CALL read_potential(qs_kind%element_symbol, potential_name, &
2356 : qs_kind%gth_potential, zeff_correction, para_env, &
2357 8199 : potential_file_name, potential_section, update_input)
2358 8199 : CALL set_potential(qs_kind%gth_potential, z=z)
2359 8199 : CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
2360 8199 : IF (.NOT. ASSOCIATED(elec_conf)) THEN
2361 8195 : CALL get_potential(potential=qs_kind%gth_potential, elec_conf=elec_conf)
2362 8195 : CALL set_qs_kind(qs_kind, elec_conf=elec_conf)
2363 : ELSE
2364 4 : CALL set_potential(potential=qs_kind%gth_potential, elec_conf=elec_conf)
2365 : END IF
2366 : CASE ("ECP")
2367 12 : CALL allocate_potential(qs_kind%sgp_potential)
2368 12 : CALL get_potential(qs_kind%sgp_potential, description=description)
2369 : CALL read_ecp_potential(ptable(z)%symbol, ecppot, &
2370 12 : potential_name, potential_file_name, potential_section)
2371 12 : IF (ecp_semi_local) THEN
2372 12 : description(1) = "Semi-local Gaussian pseudopotential "
2373 12 : description(2) = "ECP "//TRIM(potential_name)
2374 12 : description(3) = "LIBGRPP: A. V. Oleynichenko et al., Symmetry 15 197 2023"
2375 : description(4) = " "
2376 : ELSE
2377 0 : description(4) = "ECP "//TRIM(potential_name)
2378 : END IF
2379 : CALL set_potential(qs_kind%sgp_potential, name=ecppot%pname, description=description, &
2380 : zeff=ecppot%zion, z=z, ecp_local=.TRUE., ecp_semi_local=ecp_semi_local, &
2381 : nloc=ecppot%nloc, nrloc=ecppot%nrloc, aloc=ecppot%aloc, bloc=ecppot%bloc, &
2382 12 : has_nlcc=.FALSE.)
2383 : CALL set_potential(qs_kind%sgp_potential, sl_lmax=ecppot%lmax, &
2384 12 : npot=ecppot%npot, nrpot=ecppot%nrpot, apot=ecppot%apot, bpot=ecppot%bpot)
2385 : ! convert PP
2386 12 : IF (.NOT. ecp_semi_local) THEN
2387 0 : CPABORT("ECPs are only well tested in their semi-local form")
2388 0 : CALL get_qs_kind(qs_kind, basis_set=orb_basis_set)
2389 0 : CALL sgp_construction(sgp_pot=sgppot, ecp_pot=ecppot, orb_basis=orb_basis_set, error=error)
2390 0 : IF (iounit > 0) THEN
2391 0 : WRITE (iounit, "(/,T2,'PP Transformation for ',A)") TRIM(ecppot%pname)
2392 0 : IF (sgppot%has_local) THEN
2393 0 : WRITE (iounit, "(T8,'Accuracy for local part:',T41,F10.3,'%',T61,F20.12)") error(4), error(1)
2394 : END IF
2395 0 : IF (sgppot%has_nonlocal) THEN
2396 0 : WRITE (iounit, "(T8,'Accuracy for nonlocal part:',T41,F10.3,'%',T61,F20.12)") error(5), error(2)
2397 : END IF
2398 0 : IF (sgppot%has_nlcc) THEN
2399 0 : WRITE (iounit, "(T8,'Accuracy for NLCC density:',T61,F20.12)") error(3)
2400 : END IF
2401 : END IF
2402 : END IF
2403 12 : IF (sgppot%has_nonlocal) THEN
2404 : CALL set_potential(qs_kind%sgp_potential, n_nonlocal=sgppot%n_nonlocal, lmax=sgppot%lmax, &
2405 0 : is_nonlocal=sgppot%is_nonlocal)
2406 0 : nnl = sgppot%n_nonlocal
2407 0 : nppnl = 0
2408 0 : DO l = 0, sgppot%lmax
2409 0 : nppnl = nppnl + nnl*nco(l)
2410 : END DO
2411 0 : l = sgppot%lmax
2412 0 : ALLOCATE (a_nl(nnl), h_nl(nnl, 0:l), c_nl(nnl, nnl, 0:l))
2413 0 : a_nl(:) = sgppot%a_nonlocal(:)
2414 0 : h_nl(:, :) = sgppot%h_nonlocal(:, :)
2415 0 : DO l = 0, sgppot%lmax
2416 0 : c_nl(:, :, l) = sgppot%c_nonlocal(:, :, l)*SQRT(2._dp*l + 1.0_dp)
2417 : END DO
2418 0 : CALL set_potential(qs_kind%sgp_potential, nppnl=nppnl, a_nonlocal=a_nl, h_nonlocal=h_nl, c_nonlocal=c_nl)
2419 : ELSE
2420 12 : CALL set_potential(qs_kind%sgp_potential, n_nonlocal=0, lmax=-1, is_nonlocal=sgppot%is_nonlocal)
2421 12 : CALL set_potential(qs_kind%sgp_potential, nppnl=0)
2422 : END IF
2423 : !
2424 12 : CPASSERT(.NOT. sgppot%has_local)
2425 12 : CPASSERT(.NOT. sgppot%has_nlcc)
2426 : ! core
2427 12 : rc = 0.5_dp*qs_kind%covalent_radius*angstrom
2428 12 : rc = MAX(rc, 0.2_dp)
2429 12 : rc = MIN(rc, 1.0_dp)
2430 12 : alpha = 1.0_dp/(2.0_dp*rc**2)
2431 12 : ccore = ecppot%zion*SQRT((alpha/pi)**3)
2432 : CALL set_potential(qs_kind%sgp_potential, alpha_core_charge=alpha, ccore_charge=ccore, &
2433 12 : core_charge_radius=rc)
2434 12 : CALL atom_sgp_release(sgppot)
2435 12 : CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
2436 12 : IF (.NOT. ASSOCIATED(elec_conf)) THEN
2437 12 : CALL set_qs_kind(qs_kind, elec_conf=ecppot%econf)
2438 : END IF
2439 12 : CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
2440 12 : CALL set_potential(qs_kind%sgp_potential, elec_conf=elec_conf)
2441 : CASE ("UPF")
2442 12 : CALL allocate_potential(qs_kind%sgp_potential)
2443 12 : CALL get_potential(qs_kind%sgp_potential, description=description)
2444 12 : description(4) = "UPF "//TRIM(potential_name)
2445 12 : CALL atom_read_upf(upfpot, potential_name)
2446 : CALL set_potential(qs_kind%sgp_potential, name=upfpot%pname, description=description, &
2447 12 : zeff=upfpot%zion, z=z, has_nlcc=upfpot%core_correction)
2448 : ! convert pp
2449 12 : CALL sgp_construction(sgp_pot=sgppot, upf_pot=upfpot, error=error)
2450 12 : IF (iounit > 0) THEN
2451 6 : WRITE (iounit, "(/,T2,'PP Transformation for ',A)") TRIM(upfpot%pname)
2452 6 : IF (sgppot%has_local) THEN
2453 6 : WRITE (iounit, "(T8,'Accuracy for local part:',T61,F20.12)") error(1)
2454 : END IF
2455 6 : IF (sgppot%has_nonlocal) THEN
2456 3 : WRITE (iounit, "(T8,'Accuracy for nonlocal part:',T61,F20.12)") error(2)
2457 : END IF
2458 6 : IF (sgppot%has_nlcc) THEN
2459 0 : WRITE (iounit, "(T8,'Accuracy for NLCC density:',T61,F20.12)") error(3)
2460 : END IF
2461 : END IF
2462 12 : IF (sgppot%has_nonlocal) THEN
2463 : CALL set_potential(qs_kind%sgp_potential, n_nonlocal=sgppot%n_nonlocal, lmax=sgppot%lmax, &
2464 6 : is_nonlocal=sgppot%is_nonlocal)
2465 6 : nnl = sgppot%n_nonlocal
2466 6 : nppnl = 0
2467 12 : DO l = 0, sgppot%lmax
2468 12 : nppnl = nppnl + nnl*nco(l)
2469 : END DO
2470 6 : l = sgppot%lmax
2471 60 : ALLOCATE (a_nl(nnl), h_nl(nnl, 0:l), c_nl(nnl, nnl, 0:l))
2472 54 : a_nl(:) = sgppot%a_nonlocal(:)
2473 60 : h_nl(:, :) = sgppot%h_nonlocal(:, :)
2474 444 : c_nl(:, :, :) = sgppot%c_nonlocal(:, :, :)
2475 6 : CALL set_potential(qs_kind%sgp_potential, nppnl=nppnl, a_nonlocal=a_nl, h_nonlocal=h_nl, c_nonlocal=c_nl)
2476 : ELSE
2477 6 : CALL set_potential(qs_kind%sgp_potential, n_nonlocal=0, lmax=-1, is_nonlocal=sgppot%is_nonlocal)
2478 6 : CALL set_potential(qs_kind%sgp_potential, nppnl=0)
2479 : END IF
2480 12 : CPASSERT(sgppot%has_local)
2481 : ! core
2482 12 : rc = sgppot%ac_local
2483 12 : alpha = 1.0_dp/(2.0_dp*rc**2)
2484 12 : ccore = upfpot%zion*SQRT((alpha/pi)**3)
2485 : CALL set_potential(qs_kind%sgp_potential, alpha_core_charge=alpha, ccore_charge=ccore, &
2486 12 : core_charge_radius=rc)
2487 : ! local potential
2488 12 : nloc = sgppot%n_local
2489 48 : ALLOCATE (aloc(nloc), cloc(nloc))
2490 156 : aloc(1:nloc) = sgppot%a_local(1:nloc)
2491 156 : cloc(1:nloc) = sgppot%c_local(1:nloc)
2492 12 : CALL set_potential(qs_kind%sgp_potential, n_local=nloc, a_local=aloc, c_local=cloc)
2493 12 : IF (sgppot%has_nlcc) THEN
2494 0 : nlcc = sgppot%n_nlcc
2495 0 : ALLOCATE (anlcc(nlcc), cnlcc(nlcc))
2496 0 : anlcc(1:nlcc) = sgppot%a_nlcc(1:nlcc)
2497 0 : cnlcc(1:nlcc) = sgppot%c_nlcc(1:nlcc)
2498 0 : CALL set_potential(qs_kind%sgp_potential, has_nlcc=.TRUE., n_nlcc=nlcc, a_nlcc=anlcc, c_nlcc=cnlcc)
2499 : END IF
2500 12 : CALL set_potential(qs_kind%sgp_potential, z=z)
2501 12 : CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
2502 12 : IF (.NOT. ASSOCIATED(elec_conf)) THEN
2503 12 : CALL set_qs_kind(qs_kind, elec_conf=upfpot%econf)
2504 : END IF
2505 12 : CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
2506 12 : CALL set_potential(qs_kind%sgp_potential, elec_conf=elec_conf)
2507 12 : CALL atom_release_upf(upfpot)
2508 12 : CALL atom_sgp_release(sgppot)
2509 : CASE DEFAULT
2510 : CALL cp_abort(__LOCATION__, &
2511 : "An invalid potential type <"// &
2512 : TRIM(potential_name)//"> was specified "// &
2513 : "for the atomic kind <"// &
2514 9247 : TRIM(qs_kind%name))
2515 : END SELECT
2516 : ELSE
2517 : CALL cp_abort(__LOCATION__, &
2518 : "No potential type was defined for the "// &
2519 0 : "atomic kind <"//TRIM(qs_kind%name)//">")
2520 : END IF
2521 :
2522 9247 : CALL check_potential_basis_compatibility(qs_kind)
2523 :
2524 : ! Allocate and initialise the potential data set structure
2525 9247 : IF ((kgpot_name /= '') .OR. explicit_kgpot) THEN
2526 9247 : ipos = INDEX(kgpot_name, "-")
2527 9247 : IF (ipos > 1) THEN
2528 20 : kgpot_type = kgpot_name(:ipos - 1)
2529 : ELSE
2530 9227 : kgpot_type = kgpot_name
2531 : END IF
2532 9247 : CALL uppercase(kgpot_type)
2533 :
2534 9267 : SELECT CASE (TRIM(kgpot_type))
2535 : CASE ("TNADD")
2536 : ! determine the pseudopotential file to search
2537 20 : IF (kg_potential_fn_kind == "-") THEN
2538 20 : CALL section_vals_val_get(dft_section, "POTENTIAL_FILE_NAME", c_val=potential_file_name)
2539 : ELSE
2540 0 : potential_file_name = kg_potential_fn_kind
2541 : END IF
2542 20 : CALL allocate_potential(qs_kind%tnadd_potential)
2543 : CALL read_potential(qs_kind%element_symbol, kgpot_name, &
2544 : qs_kind%tnadd_potential, para_env, &
2545 20 : potential_file_name, kgpot_section, update_input)
2546 : CASE ("NONE")
2547 9227 : NULLIFY (qs_kind%tnadd_potential)
2548 : CASE DEFAULT
2549 : CALL cp_abort(__LOCATION__, &
2550 : "An invalid kg_potential type <"// &
2551 : TRIM(potential_name)//"> was specified "// &
2552 : "for the atomic kind <"// &
2553 9247 : TRIM(qs_kind%name))
2554 : END SELECT
2555 : END IF
2556 : END IF
2557 : END SELECT
2558 :
2559 12823 : CALL timestop(handle)
2560 :
2561 7668154 : END SUBROUTINE read_qs_kind
2562 :
2563 : ! **************************************************************************************************
2564 : !> \brief Ensure pseudo-potential and basis set were optimized for same number of valence electrons
2565 : !> \param qs_kind ...
2566 : !> \author Ole Schuett
2567 : ! **************************************************************************************************
2568 9247 : SUBROUTINE check_potential_basis_compatibility(qs_kind)
2569 : TYPE(qs_kind_type), INTENT(INOUT) :: qs_kind
2570 :
2571 : CHARACTER(LEN=default_string_length) :: name
2572 : INTEGER :: nbs, npp
2573 : TYPE(gth_potential_type), POINTER :: gth_potential
2574 : TYPE(gto_basis_set_type), POINTER :: basis_set
2575 :
2576 9247 : CALL get_qs_kind(qs_kind, name=name, gth_potential=gth_potential, basis_set=basis_set)
2577 :
2578 9247 : npp = -1; nbs = -1
2579 9247 : IF (ASSOCIATED(gth_potential)) &
2580 8199 : npp = parse_valence_electrons(gth_potential%aliases)
2581 9247 : IF (ASSOCIATED(basis_set)) &
2582 9247 : nbs = parse_valence_electrons(basis_set%aliases)
2583 :
2584 9247 : IF (npp >= 0 .AND. nbs >= 0 .AND. npp /= nbs) &
2585 : CALL cp_abort(__LOCATION__, "Basis-set and pseudo-potential of atomic kind '"//TRIM(name)//"'"// &
2586 0 : " were optimized for different valence electron numbers.")
2587 :
2588 9247 : END SUBROUTINE check_potential_basis_compatibility
2589 :
2590 : ! **************************************************************************************************
2591 : !> \brief Tries to parse valence eletron number using "-QXXX" notation, returns -1 if not found.
2592 : !> \param string ...
2593 : !> \return ...
2594 : !> \author Ole Schuett
2595 : ! **************************************************************************************************
2596 17446 : FUNCTION parse_valence_electrons(string) RESULT(n)
2597 : CHARACTER(*) :: string
2598 : INTEGER :: n
2599 :
2600 : INTEGER :: i, istat, j
2601 :
2602 17446 : i = INDEX(string, "-Q", .TRUE.)
2603 17446 : IF (i == 0) THEN
2604 6040 : n = -1
2605 : ELSE
2606 11406 : j = SCAN(string(i + 2:), "- ")
2607 11406 : READ (string(i + 2:i + j), '(I3)', iostat=istat) n
2608 11406 : IF (istat /= 0) n = -1
2609 : END IF
2610 :
2611 17446 : END FUNCTION
2612 :
2613 : ! **************************************************************************************************
2614 : !> \brief Read an atomic kind set data set from the input file.
2615 : !> \param qs_kind_set ...
2616 : !> \param atomic_kind_set ...
2617 : !> \param kind_section ...
2618 : !> \param para_env ...
2619 : !> \param force_env_section ...
2620 : ! **************************************************************************************************
2621 6730 : SUBROUTINE create_qs_kind_set(qs_kind_set, atomic_kind_set, kind_section, para_env, force_env_section)
2622 :
2623 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2624 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2625 : TYPE(section_vals_type), POINTER :: kind_section
2626 : TYPE(mp_para_env_type), POINTER :: para_env
2627 : TYPE(section_vals_type), POINTER :: force_env_section
2628 :
2629 : CHARACTER(len=*), PARAMETER :: routineN = 'create_qs_kind_set'
2630 :
2631 : INTEGER :: handle, ikind, method, nkind, qs_method
2632 : LOGICAL :: no_fail
2633 :
2634 6730 : CALL timeset(routineN, handle)
2635 :
2636 6730 : IF (ASSOCIATED(qs_kind_set)) CPABORT("create_qs_kind_set: qs_kind_set already associated")
2637 6730 : IF (.NOT. ASSOCIATED(atomic_kind_set)) CPABORT("create_qs_kind_set: atomic_kind_set not associated")
2638 :
2639 6730 : no_fail = .FALSE.
2640 :
2641 : ! Between all methods only SE and DFTB/xTB may not need a KIND section.
2642 6730 : CALL section_vals_val_get(force_env_section, "METHOD", i_val=method)
2643 6730 : IF (method == do_qs) THEN
2644 6710 : CALL section_vals_val_get(force_env_section, "DFT%QS%METHOD", i_val=qs_method)
2645 998 : SELECT CASE (qs_method)
2646 : CASE (do_method_mndo, do_method_am1, do_method_pm3, do_method_pm6fm, do_method_pm6, &
2647 : do_method_pdg, do_method_rm1, do_method_mndod, do_method_pnnl)
2648 998 : no_fail = .TRUE.
2649 : CASE (do_method_dftb)
2650 222 : no_fail = .TRUE.
2651 : CASE (do_method_xtb)
2652 6710 : no_fail = .TRUE.
2653 : END SELECT
2654 20 : ELSE IF (method == do_sirius) THEN
2655 16 : qs_method = do_method_pw
2656 : ELSE
2657 4 : qs_method = method
2658 : END IF
2659 :
2660 6730 : nkind = SIZE(atomic_kind_set)
2661 167613 : ALLOCATE (qs_kind_set(nkind))
2662 :
2663 19553 : DO ikind = 1, nkind
2664 12823 : qs_kind_set(ikind)%name = atomic_kind_set(ikind)%name
2665 12823 : qs_kind_set(ikind)%element_symbol = atomic_kind_set(ikind)%element_symbol
2666 12823 : qs_kind_set(ikind)%natom = atomic_kind_set(ikind)%natom
2667 19553 : CALL read_qs_kind(qs_kind_set(ikind), kind_section, para_env, force_env_section, no_fail, qs_method)
2668 : END DO
2669 :
2670 6730 : CALL timestop(handle)
2671 :
2672 13460 : END SUBROUTINE create_qs_kind_set
2673 :
2674 : ! **************************************************************************************************
2675 : !> \brief This routines should perform only checks. no settings are allowed at
2676 : !> this level anymore..
2677 : !> \param qs_kind ...
2678 : !> \param dft_control ...
2679 : !> \param subsys_section ...
2680 : ! **************************************************************************************************
2681 12745 : SUBROUTINE check_qs_kind(qs_kind, dft_control, subsys_section)
2682 :
2683 : TYPE(qs_kind_type), POINTER :: qs_kind
2684 : TYPE(dft_control_type), INTENT(IN) :: dft_control
2685 : TYPE(section_vals_type), POINTER :: subsys_section
2686 :
2687 : INTEGER :: gfn_type
2688 : LOGICAL :: defined
2689 : TYPE(qs_dftb_atom_type), POINTER :: dftb_parameter
2690 : TYPE(semi_empirical_type), POINTER :: se_parameter
2691 : TYPE(xtb_atom_type), POINTER :: xtb_parameter
2692 :
2693 12745 : IF (dft_control%qs_control%semi_empirical) THEN
2694 2240 : CALL get_qs_kind(qs_kind, se_parameter=se_parameter)
2695 2240 : CPASSERT(ASSOCIATED(se_parameter))
2696 2240 : CALL get_se_param(se_parameter, defined=defined)
2697 2240 : CPASSERT(defined)
2698 2240 : CALL write_se_param(se_parameter, subsys_section)
2699 10505 : ELSE IF (dft_control%qs_control%dftb) THEN
2700 480 : CALL get_qs_kind(qs_kind, dftb_parameter=dftb_parameter)
2701 480 : CPASSERT(ASSOCIATED(dftb_parameter))
2702 480 : CALL get_dftb_atom_param(dftb_parameter, defined=defined)
2703 480 : CPASSERT(defined)
2704 480 : CALL write_dftb_atom_param(dftb_parameter, subsys_section)
2705 10025 : ELSE IF (dft_control%qs_control%xtb) THEN
2706 692 : CALL get_qs_kind(qs_kind, xtb_parameter=xtb_parameter)
2707 692 : CPASSERT(ASSOCIATED(xtb_parameter))
2708 692 : gfn_type = dft_control%qs_control%xtb_control%gfn_type
2709 692 : CALL write_xtb_atom_param(xtb_parameter, gfn_type, subsys_section)
2710 : END IF
2711 :
2712 12745 : END SUBROUTINE check_qs_kind
2713 :
2714 : ! **************************************************************************************************
2715 : !> \brief ...
2716 : !> \param qs_kind_set ...
2717 : !> \param dft_control ...
2718 : !> \param subsys_section ...
2719 : ! **************************************************************************************************
2720 6686 : SUBROUTINE check_qs_kind_set(qs_kind_set, dft_control, subsys_section)
2721 :
2722 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2723 : TYPE(dft_control_type), INTENT(IN) :: dft_control
2724 : TYPE(section_vals_type), POINTER :: subsys_section
2725 :
2726 : CHARACTER(len=*), PARAMETER :: routineN = 'check_qs_kind_set'
2727 :
2728 : INTEGER :: handle, ikind, nkind
2729 : TYPE(qs_kind_type), POINTER :: qs_kind
2730 :
2731 6686 : CALL timeset(routineN, handle)
2732 6686 : IF (ASSOCIATED(qs_kind_set)) THEN
2733 6686 : nkind = SIZE(qs_kind_set)
2734 19431 : DO ikind = 1, nkind
2735 12745 : qs_kind => qs_kind_set(ikind)
2736 19431 : CALL check_qs_kind(qs_kind, dft_control, subsys_section)
2737 : END DO
2738 6686 : IF (dft_control%qs_control%xtb) THEN
2739 : CALL write_xtb_kab_param(qs_kind_set, subsys_section, &
2740 284 : dft_control%qs_control%xtb_control)
2741 : END IF
2742 : ELSE
2743 0 : CPABORT("The pointer qs_kind_set is not associated")
2744 : END IF
2745 6686 : CALL timestop(handle)
2746 6686 : END SUBROUTINE check_qs_kind_set
2747 :
2748 : ! **************************************************************************************************
2749 : !> \brief ...
2750 : !> \param qs_kind_set ...
2751 : !> \param subsys_section ...
2752 : !> \param xtb_control ...
2753 : ! **************************************************************************************************
2754 284 : SUBROUTINE write_xtb_kab_param(qs_kind_set, subsys_section, xtb_control)
2755 :
2756 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2757 : TYPE(section_vals_type), POINTER :: subsys_section
2758 : TYPE(xtb_control_type), POINTER :: xtb_control
2759 :
2760 : CHARACTER(LEN=default_string_length) :: aname, bname
2761 : INTEGER :: ikind, io_unit, jkind, nkind, za, zb
2762 : TYPE(cp_logger_type), POINTER :: logger
2763 : TYPE(qs_kind_type), POINTER :: qs_kinda, qs_kindb
2764 : TYPE(xtb_atom_type), POINTER :: xtb_parameter_a, xtb_parameter_b
2765 :
2766 284 : NULLIFY (logger)
2767 284 : logger => cp_get_default_logger()
2768 284 : IF (BTEST(cp_print_key_should_output(logger%iter_info, subsys_section, &
2769 : "PRINT%KINDS/POTENTIAL"), cp_p_file)) THEN
2770 :
2771 0 : io_unit = cp_print_key_unit_nr(logger, subsys_section, "PRINT%KINDS", extension=".Log")
2772 0 : IF (io_unit > 0) THEN
2773 :
2774 0 : WRITE (io_unit, "(/,T2,A)") "xTB| Kab parameters"
2775 0 : nkind = SIZE(qs_kind_set)
2776 0 : DO ikind = 1, nkind
2777 0 : qs_kinda => qs_kind_set(ikind)
2778 0 : CALL get_qs_kind(qs_kinda, xtb_parameter=xtb_parameter_a)
2779 0 : CALL get_xtb_atom_param(xtb_parameter_a, aname=aname, z=za)
2780 0 : DO jkind = ikind, nkind
2781 0 : qs_kindb => qs_kind_set(jkind)
2782 0 : CALL get_qs_kind(qs_kindb, xtb_parameter=xtb_parameter_b)
2783 0 : CALL get_xtb_atom_param(xtb_parameter_b, aname=bname, z=zb)
2784 : WRITE (io_unit, "(A,T10,A15,T25,A15,T71,F10.3)") &
2785 0 : " Kab:", TRIM(aname), TRIM(bname), xtb_set_kab(za, zb, xtb_control)
2786 : END DO
2787 : END DO
2788 0 : WRITE (io_unit, *)
2789 :
2790 : END IF
2791 :
2792 0 : CALL cp_print_key_finished_output(io_unit, logger, subsys_section, "PRINT%KINDS")
2793 : END IF
2794 :
2795 284 : END SUBROUTINE write_xtb_kab_param
2796 :
2797 : ! **************************************************************************************************
2798 : !> \brief Set the components of an atomic kind data set.
2799 : !> \param qs_kind ...
2800 : !> \param paw_atom ...
2801 : !> \param ghost ...
2802 : !> \param floating ...
2803 : !> \param hard_radius ...
2804 : !> \param hard0_radius ...
2805 : !> \param covalent_radius ...
2806 : !> \param vdw_radius ...
2807 : !> \param lmax_rho0 ...
2808 : !> \param zeff ...
2809 : !> \param no_optimize ...
2810 : !> \param dispersion ...
2811 : !> \param u_minus_j ...
2812 : !> \param reltmat ...
2813 : !> \param dftb_parameter ...
2814 : !> \param xtb_parameter ...
2815 : !> \param elec_conf ...
2816 : !> \param pao_basis_size ...
2817 : ! **************************************************************************************************
2818 20197 : SUBROUTINE set_qs_kind(qs_kind, paw_atom, ghost, floating, hard_radius, hard0_radius, &
2819 : covalent_radius, vdw_radius, lmax_rho0, zeff, &
2820 : no_optimize, dispersion, u_minus_j, reltmat, &
2821 : dftb_parameter, xtb_parameter, &
2822 20197 : elec_conf, pao_basis_size)
2823 :
2824 : TYPE(qs_kind_type), INTENT(INOUT) :: qs_kind
2825 : LOGICAL, INTENT(IN), OPTIONAL :: paw_atom, ghost, floating
2826 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: hard_radius, hard0_radius, &
2827 : covalent_radius, vdw_radius
2828 : INTEGER, INTENT(IN), OPTIONAL :: lmax_rho0
2829 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: zeff
2830 : LOGICAL, INTENT(IN), OPTIONAL :: no_optimize
2831 : TYPE(qs_atom_dispersion_type), OPTIONAL, POINTER :: dispersion
2832 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: u_minus_j
2833 : REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: reltmat
2834 : TYPE(qs_dftb_atom_type), OPTIONAL, POINTER :: dftb_parameter
2835 : TYPE(xtb_atom_type), OPTIONAL, POINTER :: xtb_parameter
2836 : INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: elec_conf
2837 : INTEGER, INTENT(IN), OPTIONAL :: pao_basis_size
2838 :
2839 20197 : IF (PRESENT(dftb_parameter)) qs_kind%dftb_parameter => dftb_parameter
2840 20197 : IF (PRESENT(xtb_parameter)) qs_kind%xtb_parameter => xtb_parameter
2841 20197 : IF (PRESENT(elec_conf)) THEN
2842 12681 : IF (ASSOCIATED(qs_kind%elec_conf)) THEN
2843 0 : DEALLOCATE (qs_kind%elec_conf)
2844 : END IF
2845 38043 : ALLOCATE (qs_kind%elec_conf(0:SIZE(elec_conf) - 1))
2846 44021 : qs_kind%elec_conf(:) = elec_conf(:)
2847 : END IF
2848 20197 : IF (PRESENT(paw_atom)) qs_kind%paw_atom = paw_atom
2849 20197 : IF (PRESENT(hard_radius)) qs_kind%hard_radius = hard_radius
2850 20197 : IF (PRESENT(hard0_radius)) qs_kind%hard0_radius = hard0_radius
2851 20197 : IF (PRESENT(covalent_radius)) qs_kind%covalent_radius = covalent_radius
2852 20197 : IF (PRESENT(vdw_radius)) qs_kind%vdw_radius = vdw_radius
2853 20197 : IF (PRESENT(lmax_rho0)) qs_kind%lmax_rho0 = lmax_rho0
2854 20197 : IF (PRESENT(zeff)) THEN
2855 0 : IF (ASSOCIATED(qs_kind%all_potential)) THEN
2856 0 : CALL set_potential(potential=qs_kind%all_potential, zeff=zeff)
2857 0 : ELSE IF (ASSOCIATED(qs_kind%gth_potential)) THEN
2858 0 : CALL set_potential(potential=qs_kind%gth_potential, zeff=zeff)
2859 0 : ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
2860 0 : CALL set_potential(potential=qs_kind%sgp_potential, zeff=zeff)
2861 : END IF
2862 : END IF
2863 20197 : IF (PRESENT(ghost)) qs_kind%ghost = ghost
2864 :
2865 20197 : IF (PRESENT(floating)) qs_kind%floating = floating
2866 :
2867 20197 : IF (PRESENT(no_optimize)) qs_kind%no_optimize = no_optimize
2868 :
2869 20197 : IF (PRESENT(dispersion)) qs_kind%dispersion => dispersion
2870 :
2871 20197 : IF (PRESENT(u_minus_j)) THEN
2872 424 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
2873 424 : qs_kind%dft_plus_u%u_minus_j = u_minus_j
2874 : END IF
2875 : END IF
2876 :
2877 20197 : IF (PRESENT(reltmat)) qs_kind%reltmat => reltmat
2878 :
2879 20197 : IF (PRESENT(pao_basis_size)) qs_kind%pao_basis_size = pao_basis_size
2880 :
2881 20197 : END SUBROUTINE set_qs_kind
2882 :
2883 : ! **************************************************************************************************
2884 : !> \brief Write an atomic kind data set to the output unit.
2885 : !> \param qs_kind ...
2886 : !> \param kind_number ...
2887 : !> \param output_unit ...
2888 : !> \par History
2889 : !> Creation (09.02.2002,MK)
2890 : ! **************************************************************************************************
2891 3556 : SUBROUTINE write_qs_kind(qs_kind, kind_number, output_unit)
2892 :
2893 : TYPE(qs_kind_type), POINTER :: qs_kind
2894 : INTEGER, INTENT(in) :: kind_number, output_unit
2895 :
2896 : CHARACTER(LEN=3) :: yon
2897 : CHARACTER(LEN=default_string_length) :: basis_type, bstring
2898 : INTEGER :: ibas
2899 : LOGICAL :: do_print
2900 : TYPE(gto_basis_set_type), POINTER :: tmp_basis
2901 :
2902 3556 : IF (output_unit > 0) THEN
2903 :
2904 3556 : IF (ASSOCIATED(qs_kind)) THEN
2905 : WRITE (UNIT=output_unit, FMT="(/,T2,I2,A,T57,A,T75,I6)") &
2906 3556 : kind_number, ". Atomic kind: "//TRIM(qs_kind%name), &
2907 7112 : "Number of atoms: ", qs_kind%natom
2908 :
2909 74676 : DO ibas = 1, SIZE(qs_kind%basis_sets, 1)
2910 71120 : NULLIFY (tmp_basis)
2911 : CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis, &
2912 71120 : inumbas=ibas, basis_type=basis_type)
2913 71120 : do_print = .TRUE.
2914 66534 : SELECT CASE (basis_type)
2915 : CASE DEFAULT
2916 66534 : bstring = "Basis Set"
2917 3470 : do_print = .FALSE.
2918 : CASE ("ORB")
2919 3470 : bstring = "Orbital Basis Set"
2920 : CASE ("ORB_SOFT")
2921 457 : bstring = "GAPW Soft Basis Set"
2922 0 : do_print = .FALSE.
2923 : CASE ("AUX")
2924 0 : bstring = "Auxiliary Basis Set"
2925 : CASE ("MIN")
2926 0 : bstring = "Minimal Basis Set"
2927 : CASE ("RI_AUX")
2928 347 : bstring = "RI Auxiliary Basis Set"
2929 : CASE ("AUX_FIT")
2930 219 : bstring = "Auxiliary Fit Basis Set"
2931 : CASE ("LRI_AUX")
2932 15 : bstring = "LRI Basis Set"
2933 : CASE ("P_LRI_AUX")
2934 4 : bstring = "LRI Basis Set for TDDFPT"
2935 : CASE ("RI_XAS")
2936 0 : bstring = "RI XAS Basis Set"
2937 : CASE ("RI_HFX")
2938 71120 : bstring = "RI HFX Basis Set"
2939 : END SELECT
2940 :
2941 3556 : IF (do_print) THEN
2942 4129 : CALL write_orb_basis_set(tmp_basis, output_unit, bstring)
2943 : END IF
2944 :
2945 : END DO
2946 :
2947 3556 : IF (qs_kind%ghost) THEN
2948 : WRITE (UNIT=output_unit, FMT="(/,T6,A)") &
2949 7 : "The atoms of this atomic kind are GHOST atoms!"
2950 : END IF
2951 3556 : IF (qs_kind%floating) THEN
2952 : WRITE (UNIT=output_unit, FMT="(/,T6,A)") &
2953 0 : "The atoms of this atomic kind are FLOATING BASIS FUNCTIONS."
2954 : END IF
2955 3556 : IF (qs_kind%covalent_radius > 0.0_dp) THEN
2956 : WRITE (UNIT=output_unit, FMT="(/,T8,A,T71,F10.3)") &
2957 2396 : "Atomic covalent radius [Angstrom]:", &
2958 4792 : qs_kind%covalent_radius*angstrom
2959 : END IF
2960 3556 : IF (qs_kind%vdw_radius > 0.0_dp) THEN
2961 : WRITE (UNIT=output_unit, FMT="(/,T8,A,T71,F10.3)") &
2962 2396 : "Atomic van der Waals radius [Angstrom]:", &
2963 4792 : qs_kind%vdw_radius*angstrom
2964 : END IF
2965 3556 : IF (qs_kind%paw_atom) THEN
2966 : WRITE (UNIT=output_unit, FMT="(/,T6,A)") &
2967 368 : "The atoms of this atomic kind are PAW atoms (GAPW):"
2968 : WRITE (UNIT=output_unit, FMT="(T8,A,T71,F10.3)") &
2969 368 : "Hard Gaussian function radius:", qs_kind%hard_radius, &
2970 368 : "Rho0 radius:", qs_kind%hard0_radius, &
2971 368 : "Maximum GTO radius used for PAW projector construction:", &
2972 736 : qs_kind%max_rad_local
2973 368 : NULLIFY (tmp_basis)
2974 : CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis, &
2975 368 : basis_type="ORB_SOFT")
2976 368 : CALL write_orb_basis_set(tmp_basis, output_unit, "GAPW Soft Basis Set")
2977 : END IF
2978 : ! Potentials
2979 3556 : IF (ASSOCIATED(qs_kind%all_potential)) CALL write_potential(qs_kind%all_potential, output_unit)
2980 3556 : IF (ASSOCIATED(qs_kind%gth_potential)) CALL write_potential(qs_kind%gth_potential, output_unit)
2981 3556 : IF (ASSOCIATED(qs_kind%sgp_potential)) CALL write_potential(qs_kind%sgp_potential, output_unit)
2982 3556 : IF (ASSOCIATED(qs_kind%tnadd_potential)) CALL write_potential(qs_kind%tnadd_potential, output_unit)
2983 3556 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
2984 : WRITE (UNIT=output_unit, FMT="(/,T6,A,/,T8,A,T76,I5,/,T8,A,T73,F8.3)") &
2985 16 : "A DFT+U correction is applied to atoms of this atomic kind:", &
2986 16 : "Angular quantum momentum number L:", qs_kind%dft_plus_u%l, &
2987 32 : "U(eff) = (U - J) value in [eV]:", qs_kind%dft_plus_u%u_minus_j_target*evolt
2988 16 : IF (qs_kind%dft_plus_u%u_ramping > 0.0_dp) THEN
2989 4 : IF (qs_kind%dft_plus_u%init_u_ramping_each_scf) THEN
2990 2 : yon = "YES"
2991 : ELSE
2992 2 : yon = " NO"
2993 : END IF
2994 : WRITE (UNIT=output_unit, FMT="(T8,A,T73,F8.3,/,T8,A,T73,ES8.1,/,T8,A,T78,A3)") &
2995 4 : "Increment for U ramping in [eV]:", qs_kind%dft_plus_u%u_ramping*evolt, &
2996 4 : "SCF threshold value for U ramping:", qs_kind%dft_plus_u%eps_u_ramping, &
2997 8 : "Set U ramping value to zero before each wavefunction optimisation:", yon
2998 : END IF
2999 16 : IF (ASSOCIATED(qs_kind%dft_plus_u%orbitals)) THEN
3000 : WRITE (UNIT=output_unit, FMT="(T8,A)") &
3001 2 : "An initial orbital occupation is requested:"
3002 2 : IF (ASSOCIATED(qs_kind%dft_plus_u%nelec)) THEN
3003 4 : IF (ANY(qs_kind%dft_plus_u%nelec(:) >= 0.5_dp)) THEN
3004 0 : IF (SIZE(qs_kind%dft_plus_u%nelec) > 1) THEN
3005 : WRITE (UNIT=output_unit, FMT="(T9,A,T75,F6.2)") &
3006 0 : "Number of alpha electrons:", &
3007 0 : qs_kind%dft_plus_u%nelec(1), &
3008 0 : "Number of beta electrons:", &
3009 0 : qs_kind%dft_plus_u%nelec(2)
3010 : ELSE
3011 : WRITE (UNIT=output_unit, FMT="(T9,A,T75,F6.2)") &
3012 0 : "Number of electrons:", &
3013 0 : qs_kind%dft_plus_u%nelec(1)
3014 : END IF
3015 : END IF
3016 : END IF
3017 : WRITE (UNIT=output_unit, FMT="(T9,A,(T78,I3))") &
3018 2 : "Preferred (initial) orbital occupation order (orbital M values):", &
3019 4 : qs_kind%dft_plus_u%orbitals(:)
3020 : WRITE (UNIT=output_unit, FMT="(T9,A,T71,ES10.3,/,T9,A,T76,I5)") &
3021 2 : "Threshold value for the SCF convergence criterion:", &
3022 2 : qs_kind%dft_plus_u%eps_scf, &
3023 2 : "Number of initial SCF iterations:", &
3024 4 : qs_kind%dft_plus_u%max_scf
3025 2 : IF (qs_kind%dft_plus_u%smear) THEN
3026 : WRITE (UNIT=output_unit, FMT="(T9,A)") &
3027 2 : "A smearing of the orbital occupations will be performed"
3028 : END IF
3029 : END IF
3030 : END IF
3031 : ELSE
3032 0 : CPABORT("")
3033 : END IF
3034 :
3035 : END IF
3036 :
3037 3556 : END SUBROUTINE write_qs_kind
3038 :
3039 : ! **************************************************************************************************
3040 : !> \brief Write an atomic kind set data set to the output unit.
3041 : !> \param qs_kind_set ...
3042 : !> \param subsys_section ...
3043 : !> \par History
3044 : !> Creation (09.02.2002,MK)
3045 : ! **************************************************************************************************
3046 6700 : SUBROUTINE write_qs_kind_set(qs_kind_set, subsys_section)
3047 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3048 : TYPE(section_vals_type), POINTER :: subsys_section
3049 :
3050 : CHARACTER(len=*), PARAMETER :: routineN = 'write_qs_kind_set'
3051 :
3052 : INTEGER :: handle, ikind, nkind, output_unit
3053 : TYPE(cp_logger_type), POINTER :: logger
3054 : TYPE(qs_kind_type), POINTER :: qs_kind
3055 :
3056 6700 : CALL timeset(routineN, handle)
3057 :
3058 6700 : NULLIFY (logger)
3059 6700 : logger => cp_get_default_logger()
3060 : output_unit = cp_print_key_unit_nr(logger, subsys_section, &
3061 6700 : "PRINT%KINDS", extension=".Log")
3062 6700 : IF (output_unit > 0) THEN
3063 1898 : IF (ASSOCIATED(qs_kind_set)) THEN
3064 1898 : WRITE (UNIT=output_unit, FMT="(/,/,T2,A)") "ATOMIC KIND INFORMATION"
3065 1898 : nkind = SIZE(qs_kind_set)
3066 5454 : DO ikind = 1, nkind
3067 3556 : qs_kind => qs_kind_set(ikind)
3068 5454 : CALL write_qs_kind(qs_kind, ikind, output_unit)
3069 : END DO
3070 : ELSE
3071 0 : CPABORT("")
3072 : END IF
3073 : END IF
3074 :
3075 : CALL cp_print_key_finished_output(output_unit, logger, subsys_section, &
3076 6700 : "PRINT%KINDS")
3077 :
3078 6700 : CALL timestop(handle)
3079 :
3080 6700 : END SUBROUTINE write_qs_kind_set
3081 :
3082 : ! **************************************************************************************************
3083 : !> \brief Write all the GTO basis sets of an atomic kind set to the output
3084 : !> unit (for the printing of the unnormalized basis sets as read from
3085 : !> database).
3086 : !> \param qs_kind_set ...
3087 : !> \param subsys_section ...
3088 : !> \par History
3089 : !> Creation (17.01.2002,MK)
3090 : ! **************************************************************************************************
3091 6684 : SUBROUTINE write_gto_basis_sets(qs_kind_set, subsys_section)
3092 :
3093 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3094 : TYPE(section_vals_type), POINTER :: subsys_section
3095 :
3096 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_gto_basis_sets'
3097 :
3098 : CHARACTER(LEN=default_string_length) :: basis_type, bstring
3099 : INTEGER :: handle, ibas, ikind, nkind, output_unit
3100 : TYPE(cp_logger_type), POINTER :: logger
3101 : TYPE(gto_basis_set_type), POINTER :: tmp_basis
3102 : TYPE(qs_kind_type), POINTER :: qs_kind
3103 :
3104 6684 : CALL timeset(routineN, handle)
3105 :
3106 6684 : NULLIFY (logger)
3107 6684 : logger => cp_get_default_logger()
3108 : output_unit = cp_print_key_unit_nr(logger, subsys_section, &
3109 : "PRINT%KINDS/BASIS_SET", &
3110 6684 : extension=".Log")
3111 6684 : IF (output_unit > 0) THEN
3112 60 : IF (ASSOCIATED(qs_kind_set)) THEN
3113 : WRITE (UNIT=output_unit, FMT="(/,/,T2,A)") &
3114 60 : "BASIS SET INFORMATION (Unnormalised Gaussian-type functions)"
3115 60 : nkind = SIZE(qs_kind_set)
3116 175 : DO ikind = 1, nkind
3117 115 : qs_kind => qs_kind_set(ikind)
3118 : WRITE (UNIT=output_unit, FMT="(/,T2,I2,A)") &
3119 115 : ikind, ". Atomic kind: "//TRIM(qs_kind%name)
3120 :
3121 2475 : DO ibas = 1, SIZE(qs_kind%basis_sets, 1)
3122 2300 : NULLIFY (tmp_basis)
3123 : CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis, &
3124 2300 : inumbas=ibas, basis_type=basis_type)
3125 2300 : IF (basis_type == "") CYCLE
3126 11 : SELECT CASE (basis_type)
3127 : CASE DEFAULT
3128 11 : bstring = "Basis Set"
3129 : CASE ("ORB")
3130 115 : bstring = "Orbital Basis Set"
3131 : CASE ("ORB_SOFT")
3132 11 : bstring = "GAPW Soft Basis Set"
3133 : CASE ("AUX")
3134 0 : bstring = "Auxiliary Basis Set"
3135 : CASE ("MIN")
3136 0 : bstring = "Minimal Basis Set"
3137 : CASE ("RI_AUX")
3138 0 : bstring = "RI Auxiliary Basis Set"
3139 : CASE ("AUX_FIT")
3140 0 : bstring = "Auxiliary Fit Basis Set"
3141 : CASE ("LRI_AUX")
3142 2 : bstring = "LRI Basis Set"
3143 : CASE ("P_LRI_AUX")
3144 0 : bstring = "LRI Basis Set for TDDFPT"
3145 : CASE ("RI_HFX")
3146 139 : bstring = "RI HFX Basis Set"
3147 : END SELECT
3148 :
3149 254 : IF (ASSOCIATED(tmp_basis)) CALL write_gto_basis_set(tmp_basis, output_unit, bstring)
3150 :
3151 : END DO
3152 :
3153 : END DO
3154 : ELSE
3155 0 : CPABORT("")
3156 : END IF
3157 : END IF
3158 :
3159 : CALL cp_print_key_finished_output(output_unit, logger, subsys_section, &
3160 6684 : "PRINT%KINDS/BASIS_SET")
3161 :
3162 6684 : CALL timestop(handle)
3163 :
3164 6684 : END SUBROUTINE write_gto_basis_sets
3165 :
3166 : ! **************************************************************************************************
3167 : !> \brief ...
3168 : !> \param atomic_kind ...
3169 : !> \param qs_kind ...
3170 : !> \param ncalc ...
3171 : !> \param ncore ...
3172 : !> \param nelem ...
3173 : !> \param edelta ...
3174 : ! **************************************************************************************************
3175 84406 : SUBROUTINE init_atom_electronic_state(atomic_kind, qs_kind, ncalc, ncore, nelem, edelta)
3176 :
3177 : TYPE(atomic_kind_type), INTENT(IN) :: atomic_kind
3178 : TYPE(qs_kind_type), INTENT(IN) :: qs_kind
3179 : INTEGER, DIMENSION(0:lmat, 10), INTENT(OUT) :: ncalc, ncore, nelem
3180 : REAL(KIND=dp), DIMENSION(0:lmat, 10, 2), &
3181 : INTENT(OUT) :: edelta
3182 :
3183 : INTEGER :: i, ii, is, l, ll, ne, nn, z
3184 42203 : INTEGER, DIMENSION(:), POINTER :: econf
3185 42203 : INTEGER, DIMENSION(:, :), POINTER :: addel, laddel, naddel
3186 : LOGICAL :: bs_occupation
3187 : REAL(KIND=dp) :: dmag, magnetization
3188 : TYPE(gth_potential_type), POINTER :: gth_potential
3189 : TYPE(sgp_potential_type), POINTER :: sgp_potential
3190 :
3191 42203 : CALL get_atomic_kind(atomic_kind, z=z)
3192 42203 : NULLIFY (gth_potential)
3193 : CALL get_qs_kind(qs_kind, &
3194 : gth_potential=gth_potential, &
3195 : sgp_potential=sgp_potential, &
3196 : magnetization=magnetization, &
3197 : bs_occupation=bs_occupation, &
3198 42203 : addel=addel, laddel=laddel, naddel=naddel)
3199 :
3200 : ! electronic state
3201 42203 : nelem = 0
3202 42203 : ncore = 0
3203 42203 : ncalc = 0
3204 42203 : edelta = 0.0_dp
3205 42203 : IF (ASSOCIATED(gth_potential)) THEN
3206 24061 : CALL get_potential(gth_potential, elec_conf=econf)
3207 24061 : CALL set_pseudo_state(econf, z, ncalc, ncore, nelem)
3208 18142 : ELSE IF (ASSOCIATED(sgp_potential)) THEN
3209 58 : CALL get_potential(sgp_potential, elec_conf=econf)
3210 58 : CALL set_pseudo_state(econf, z, ncalc, ncore, nelem)
3211 : ELSE
3212 90420 : DO l = 0, MIN(lmat, UBOUND(ptable(z)%e_conv, 1))
3213 72336 : ll = 2*(2*l + 1)
3214 72336 : nn = ptable(z)%e_conv(l)
3215 72336 : ii = 0
3216 18084 : DO
3217 102568 : ii = ii + 1
3218 102568 : IF (nn <= ll) THEN
3219 72336 : nelem(l, ii) = nn
3220 : EXIT
3221 : ELSE
3222 30232 : nelem(l, ii) = ll
3223 30232 : nn = nn - ll
3224 : END IF
3225 : END DO
3226 : END DO
3227 1283964 : ncalc = nelem - ncore
3228 : END IF
3229 :
3230 : ! readjust the occupation number of the orbitals as requested by user
3231 : ! this is done to break symmetry (bs) and bias the initial guess
3232 : ! to the pre-defined multiplicity/charge state of the atom
3233 42203 : IF (bs_occupation) THEN
3234 636 : DO is = 1, 2
3235 1156 : DO i = 1, SIZE(addel, 1)
3236 520 : ne = addel(i, is)
3237 520 : l = laddel(i, is)
3238 520 : nn = naddel(i, is) - l
3239 944 : IF (ne /= 0) THEN
3240 492 : IF (nn == 0) THEN
3241 0 : DO ii = SIZE(nelem, 2), 1, -1
3242 0 : IF (ncalc(l, ii) > 0) THEN
3243 0 : IF ((ncalc(l, ii) + ne) < 2*(2*l + 1) + 1) THEN
3244 0 : edelta(l, ii, is) = edelta(l, ii, is) + ne
3245 0 : nn = ii
3246 : ELSE
3247 0 : edelta(l, ii + 1, is) = edelta(l, ii + 1, is) + ne
3248 0 : nn = ii + 1
3249 : END IF
3250 : EXIT
3251 0 : ELSE IF (ii == 1) THEN
3252 0 : edelta(l, ii, is) = edelta(l, ii, is) + ne
3253 0 : nn = ii
3254 : END IF
3255 : END DO
3256 : ELSE
3257 492 : edelta(l, nn, is) = edelta(l, nn, is) + ne
3258 : END IF
3259 492 : IF (ncalc(l, nn) + edelta(l, nn, is) < 0) THEN
3260 0 : edelta(l, nn, is) = -ncalc(l, nn)
3261 : END IF
3262 : END IF
3263 : END DO
3264 : END DO
3265 30316 : edelta = 0.5_dp*edelta
3266 41991 : ELSE IF (magnetization /= 0.0_dp) THEN
3267 0 : dmag = 0.5_dp*ABS(magnetization)
3268 0 : DO l = 0, MIN(lmat, UBOUND(ptable(z)%e_conv, 1))
3269 0 : ll = 2*(2*l + 1)
3270 0 : ii = 0
3271 0 : DO i = 1, SIZE(ncalc, 2)
3272 0 : IF (ncalc(l, i) == 0) CYCLE
3273 0 : IF (ncalc(l, i) == ll) CYCLE
3274 0 : IF (ncalc(l, i) > dmag .AND. (ll - ncalc(l, i)) > dmag) THEN
3275 : ii = i
3276 : EXIT
3277 : END IF
3278 : END DO
3279 0 : IF (ii /= 0) THEN
3280 0 : edelta(l, ii, 1) = magnetization*0.5_dp
3281 0 : edelta(l, ii, 2) = -magnetization*0.5_dp
3282 0 : EXIT
3283 : END IF
3284 : END DO
3285 0 : IF (ii == 0) THEN
3286 : CALL cp_abort(__LOCATION__, &
3287 0 : "Magnetization value cannot be imposed for this atom type")
3288 : END IF
3289 : END IF
3290 :
3291 42203 : IF (qs_kind%ghost .OR. qs_kind%floating) THEN
3292 384 : nelem = 0
3293 384 : ncore = 0
3294 384 : ncalc = 0
3295 384 : edelta = 0.0_dp
3296 : END IF
3297 :
3298 42203 : END SUBROUTINE init_atom_electronic_state
3299 :
3300 : ! **************************************************************************************************
3301 : !> \brief ...
3302 : !> \param econf ...
3303 : !> \param z ...
3304 : !> \param ncalc ...
3305 : !> \param ncore ...
3306 : !> \param nelem ...
3307 : ! **************************************************************************************************
3308 24173 : SUBROUTINE set_pseudo_state(econf, z, ncalc, ncore, nelem)
3309 : INTEGER, DIMENSION(:), POINTER :: econf
3310 : INTEGER, INTENT(IN) :: z
3311 : INTEGER, DIMENSION(0:lmat, 10), INTENT(OUT) :: ncalc, ncore, nelem
3312 :
3313 : CHARACTER(LEN=default_string_length) :: message
3314 : INTEGER :: ii, iounit, l, ll, lmin, nc, nn
3315 : INTEGER, DIMENSION(0:lmat) :: econfx
3316 : TYPE(cp_logger_type), POINTER :: logger
3317 :
3318 24173 : NULLIFY (logger)
3319 24173 : logger => cp_get_default_logger()
3320 24173 : iounit = cp_logger_get_default_io_unit(logger)
3321 :
3322 24173 : econfx = 0
3323 66196 : econfx(0:SIZE(econf) - 1) = econf
3324 66196 : IF (SUM(econf) >= 0) THEN
3325 66128 : lmin = MIN(lmat, UBOUND(ptable(z)%e_conv, 1))
3326 : ! number of core electrons
3327 66128 : nc = z - SUM(econf)
3328 : ! setup ncore
3329 24139 : ncore = 0
3330 9014 : SELECT CASE (nc)
3331 : CASE (0)
3332 : CASE (2)
3333 9014 : ncore(0, 1) = 2
3334 : CASE (10)
3335 2232 : ncore(0, 1) = 2
3336 2232 : ncore(0, 2) = 2
3337 2232 : ncore(1, 1) = 6
3338 : CASE (18)
3339 58 : ncore(0, 1) = 2
3340 58 : ncore(0, 2) = 2
3341 58 : ncore(0, 3) = 2
3342 58 : ncore(1, 1) = 6
3343 58 : ncore(1, 2) = 6
3344 : CASE (28)
3345 8 : ncore(0, 1) = 2
3346 8 : ncore(0, 2) = 2
3347 8 : ncore(0, 3) = 2
3348 8 : ncore(1, 1) = 6
3349 8 : ncore(1, 2) = 6
3350 8 : ncore(2, 1) = 10
3351 : CASE (36)
3352 0 : ncore(0, 1) = 2
3353 0 : ncore(0, 2) = 2
3354 0 : ncore(0, 3) = 2
3355 0 : ncore(0, 4) = 2
3356 0 : ncore(1, 1) = 6
3357 0 : ncore(1, 2) = 6
3358 0 : ncore(1, 3) = 6
3359 0 : ncore(2, 1) = 10
3360 : CASE (46)
3361 48 : ncore(0, 1) = 2
3362 48 : ncore(0, 2) = 2
3363 48 : ncore(0, 3) = 2
3364 48 : ncore(0, 4) = 2
3365 48 : ncore(1, 1) = 6
3366 48 : ncore(1, 2) = 6
3367 48 : ncore(1, 3) = 6
3368 48 : ncore(2, 1) = 10
3369 48 : ncore(2, 2) = 10
3370 : CASE (54)
3371 4 : ncore(0, 1) = 2
3372 4 : ncore(0, 2) = 2
3373 4 : ncore(0, 3) = 2
3374 4 : ncore(0, 4) = 2
3375 4 : ncore(0, 5) = 2
3376 4 : ncore(1, 1) = 6
3377 4 : ncore(1, 2) = 6
3378 4 : ncore(1, 3) = 6
3379 4 : ncore(1, 4) = 6
3380 4 : ncore(2, 1) = 10
3381 4 : ncore(2, 2) = 10
3382 : CASE (60)
3383 18 : ncore(0, 1) = 2
3384 18 : ncore(0, 2) = 2
3385 18 : ncore(0, 3) = 2
3386 18 : ncore(0, 4) = 2
3387 18 : ncore(1, 1) = 6
3388 18 : ncore(1, 2) = 6
3389 18 : ncore(1, 3) = 6
3390 18 : ncore(2, 1) = 10
3391 18 : ncore(2, 2) = 10
3392 18 : ncore(3, 1) = 14
3393 : CASE (68)
3394 172 : ncore(0, 1) = 2
3395 172 : ncore(0, 2) = 2
3396 172 : ncore(0, 3) = 2
3397 172 : ncore(0, 4) = 2
3398 172 : ncore(0, 5) = 2
3399 172 : ncore(1, 1) = 6
3400 172 : ncore(1, 2) = 6
3401 172 : ncore(1, 3) = 6
3402 172 : ncore(1, 4) = 6
3403 172 : ncore(2, 1) = 10
3404 172 : ncore(2, 2) = 10
3405 172 : ncore(3, 1) = 14
3406 : CASE (78)
3407 12 : ncore(0, 1) = 2
3408 12 : ncore(0, 2) = 2
3409 12 : ncore(0, 3) = 2
3410 12 : ncore(0, 4) = 2
3411 12 : ncore(0, 5) = 2
3412 12 : ncore(1, 1) = 6
3413 12 : ncore(1, 2) = 6
3414 12 : ncore(1, 3) = 6
3415 12 : ncore(1, 4) = 6
3416 12 : ncore(2, 1) = 10
3417 12 : ncore(2, 2) = 10
3418 12 : ncore(2, 3) = 10
3419 12 : ncore(3, 1) = 14
3420 : CASE DEFAULT
3421 24139 : ncore(0, 1) = -1
3422 : END SELECT
3423 : ! special cases of double assignments
3424 24139 : IF (z == 65 .AND. econfx(3) == 0) THEN
3425 : ! 4f in core for Tb
3426 4 : ncore = 0
3427 4 : ncore(0, 1) = -1
3428 : END IF
3429 : ! if there is still no core, check for special cases
3430 24139 : IF (ncore(0, 1) <= 0) THEN
3431 12577 : IF (z >= 58 .AND. z <= 71) THEN
3432 : ! 4f-in-core PPs for lanthanides
3433 280 : nc = z - SUM(econf)
3434 : ! setup ncore
3435 56 : ncore = 0
3436 0 : SELECT CASE (nc)
3437 : CASE (29:42)
3438 0 : ncore(0, 1) = 2
3439 0 : ncore(0, 2) = 2
3440 0 : ncore(0, 3) = 2
3441 0 : ncore(1, 1) = 6
3442 0 : ncore(1, 2) = 6
3443 0 : ncore(2, 1) = 10
3444 0 : ncore(3, 1) = nc - 28
3445 : message = "A small-core pseudopotential with 4f-in-core is used for the lanthanide "// &
3446 0 : TRIM(ptable(z)%symbol)
3447 0 : CPHINT(TRIM(message))
3448 : CASE (47:60)
3449 56 : ncore(0, 1) = 2
3450 56 : ncore(0, 2) = 2
3451 56 : ncore(0, 3) = 2
3452 56 : ncore(0, 4) = 2
3453 56 : ncore(1, 1) = 6
3454 56 : ncore(1, 2) = 6
3455 56 : ncore(1, 3) = 6
3456 56 : ncore(2, 1) = 10
3457 56 : ncore(2, 2) = 10
3458 56 : ncore(3, 1) = nc - 46
3459 : message = "A medium-core pseudopotential with 4f-in-core is used for the lanthanide "// &
3460 56 : TRIM(ptable(z)%symbol)
3461 56 : CPHINT(TRIM(message))
3462 : CASE DEFAULT
3463 56 : ncore(0, 1) = -1
3464 : END SELECT
3465 : END IF
3466 : END IF
3467 : ! if the core is established, finish the setup
3468 24139 : IF (ncore(0, 1) >= 0) THEN
3469 120695 : DO l = 0, lmin
3470 96556 : ll = 2*(2*l + 1)
3471 1062116 : nn = SUM(ncore(l, :)) + econfx(l)
3472 96556 : ii = 0
3473 24139 : DO
3474 115672 : ii = ii + 1
3475 115672 : IF (nn <= ll) THEN
3476 96556 : nelem(l, ii) = nn
3477 : EXIT
3478 : ELSE
3479 19116 : nelem(l, ii) = ll
3480 19116 : nn = nn - ll
3481 : END IF
3482 : END DO
3483 : END DO
3484 1713869 : ncalc = nelem - ncore
3485 : ELSE
3486 : ! test for compatibility of valence occupation and full atomic occupation
3487 0 : IF (iounit > 0) THEN
3488 0 : WRITE (iounit, "(/,A,A2)") "WARNING: Core states irregular for atom type ", ptable(z)%symbol
3489 0 : WRITE (iounit, "(A,10I3)") "WARNING: Redefine ELEC_CONF in the KIND section"
3490 0 : CPABORT("Incompatible Atomic Occupations Detected")
3491 : END IF
3492 : END IF
3493 : ELSE
3494 34 : lmin = MIN(lmat, UBOUND(ptable(z)%e_conv, 1))
3495 34 : ncore = 0
3496 34 : ncalc = 0
3497 170 : DO l = 0, lmin
3498 136 : ll = 2*(2*l + 1)
3499 136 : nn = ABS(econfx(l))
3500 136 : ii = 0
3501 34 : DO
3502 136 : ii = ii + 1
3503 136 : IF (nn <= ll) THEN
3504 136 : ncalc(l, ii) = -nn
3505 : EXIT
3506 : ELSE
3507 0 : ncalc(l, ii) = -ll
3508 0 : nn = nn - ll
3509 : END IF
3510 : END DO
3511 : END DO
3512 34 : nelem = ncalc
3513 : END IF
3514 :
3515 24173 : END SUBROUTINE set_pseudo_state
3516 :
3517 : ! **************************************************************************************************
3518 : !> \brief finds if a given qs run needs to use nlcc
3519 : !> \param qs_kind_set ...
3520 : !> \return ...
3521 : ! **************************************************************************************************
3522 27044 : FUNCTION has_nlcc(qs_kind_set) RESULT(nlcc)
3523 :
3524 : TYPE(qs_kind_type), DIMENSION(:) :: qs_kind_set
3525 : LOGICAL :: nlcc
3526 :
3527 : INTEGER :: ikind
3528 : LOGICAL :: nlcc_present
3529 : TYPE(gth_potential_type), POINTER :: gth_potential
3530 : TYPE(sgp_potential_type), POINTER :: sgp_potential
3531 :
3532 27044 : nlcc = .FALSE.
3533 :
3534 80721 : DO ikind = 1, SIZE(qs_kind_set)
3535 53677 : CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential, sgp_potential=sgp_potential)
3536 80721 : IF (ASSOCIATED(gth_potential)) THEN
3537 34635 : CALL get_potential(potential=gth_potential, nlcc_present=nlcc_present)
3538 34635 : nlcc = nlcc .OR. nlcc_present
3539 19042 : ELSEIF (ASSOCIATED(sgp_potential)) THEN
3540 272 : CALL get_potential(potential=sgp_potential, has_nlcc=nlcc_present)
3541 272 : nlcc = nlcc .OR. nlcc_present
3542 : END IF
3543 : END DO
3544 :
3545 27044 : END FUNCTION has_nlcc
3546 :
3547 : ! **************************************************************************************************
3548 :
3549 0 : END MODULE qs_kind_types
|
