Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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, plus_u_lowdin, plus_u_tensorial
65 : USE input_section_types, ONLY: section_vals_get,&
66 : section_vals_get_subs_vals,&
67 : section_vals_get_subs_vals2,&
68 : section_vals_type,&
69 : section_vals_val_get
70 : USE kinds, ONLY: default_path_length,&
71 : default_string_length,&
72 : dp
73 : USE mathconstants, ONLY: pi
74 : USE message_passing, ONLY: mp_para_env_type
75 : USE orbital_pointers, ONLY: init_orbital_pointers,&
76 : nco,&
77 : ncoset
78 : USE paw_proj_set_types, ONLY: allocate_paw_proj_set,&
79 : deallocate_paw_proj_set,&
80 : get_paw_proj_set,&
81 : paw_proj_set_type,&
82 : projectors
83 : USE periodic_table, ONLY: get_ptable_info,&
84 : ptable
85 : USE physcon, ONLY: angstrom,&
86 : bohr,&
87 : evolt
88 : USE qs_cneo_types, ONLY: allocate_cneo_potential,&
89 : cneo_potential_type,&
90 : deallocate_cneo_potential,&
91 : get_cneo_potential,&
92 : set_cneo_potential,&
93 : write_cneo_potential
94 : USE qs_dftb_types, ONLY: qs_dftb_atom_type
95 : USE qs_dftb_utils, ONLY: deallocate_dftb_atom_param,&
96 : get_dftb_atom_param,&
97 : write_dftb_atom_param
98 : USE qs_dispersion_types, ONLY: qs_atom_dispersion_type
99 : USE qs_grid_atom, ONLY: allocate_grid_atom,&
100 : deallocate_grid_atom,&
101 : grid_atom_type
102 : USE qs_harmonics_atom, ONLY: allocate_harmonics_atom,&
103 : deallocate_harmonics_atom,&
104 : harmonics_atom_type
105 : USE semi_empirical_types, ONLY: get_se_param,&
106 : semi_empirical_create,&
107 : semi_empirical_release,&
108 : semi_empirical_type,&
109 : write_se_param
110 : USE semi_empirical_utils, ONLY: init_se_param,&
111 : se_param_set_default
112 : USE soft_basis_set, ONLY: create_soft_basis
113 : USE string_utilities, ONLY: uppercase
114 : USE xc_input_constants, ONLY: skala_gapw_cp2k_default,&
115 : skala_gapw_direct_valence,&
116 : skala_gapw_paw_one_center,&
117 : skala_gapw_paw_one_center_split
118 : USE xtb_parameters, ONLY: xtb_set_kab
119 : USE xtb_types, ONLY: deallocate_xtb_atom_param,&
120 : get_xtb_atom_param,&
121 : write_xtb_atom_param,&
122 : xtb_atom_type
123 : #include "./base/base_uses.f90"
124 :
125 : IMPLICIT NONE
126 :
127 : PRIVATE
128 :
129 : ! Global parameters (only in this module)
130 :
131 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_kind_types'
132 :
133 : ! **************************************************************************************************
134 : !> \brief Input parameters for the DFT+U method
135 : ! **************************************************************************************************
136 : TYPE dft_plus_u_type
137 : INTEGER :: l = -1
138 : INTEGER :: n = -1
139 : INTEGER :: max_scf = -1
140 : INTEGER :: lr_atom = -1
141 : REAL(KIND=dp) :: eps_u_ramping = 0.0_dp
142 : REAL(KIND=dp) :: eps_scf = HUGE(0.0_dp)
143 : REAL(KIND=dp) :: u_minus_j_target = 0.0_dp
144 : REAL(KIND=dp) :: u_minus_j = 0.0_dp
145 : REAL(KIND=dp) :: hund_j = 0.0_dp
146 : REAL(KIND=dp) :: u_ramping = 0.0_dp
147 : REAL(KIND=dp) :: U = 0.0_dp
148 : REAL(KIND=dp) :: J = 0.0_dp
149 : REAL(KIND=dp) :: alpha = 0.0_dp
150 : REAL(KIND=dp) :: beta = 0.0_dp
151 : REAL(KIND=dp) :: J0 = 0.0_dp
152 : REAL(KIND=dp) :: occupation = -1.0_dp
153 : REAL(KIND=dp) :: proj_shell_charge(4) = 0.0_dp
154 : REAL(KIND=dp), DIMENSION(:), &
155 : ALLOCATABLE :: perturbation_strength
156 : REAL(KIND=dp), DIMENSION(:), POINTER :: ao_coef => Null()
157 : INTEGER, DIMENSION(:), POINTER :: orbitals => Null()
158 : LOGICAL :: init_u_ramping_each_scf = .FALSE.
159 : LOGICAL :: smear = .FALSE.
160 : LOGICAL :: do_mtlr = .FALSE.
161 : LOGICAL :: u_j_loop = .FALSE.
162 : REAL(KIND=dp), DIMENSION(:), POINTER :: nelec => Null()
163 : END TYPE dft_plus_u_type
164 :
165 : ! **************************************************************************************************
166 : !> \brief Holds information about a PAO potential
167 : ! **************************************************************************************************
168 : TYPE pao_potential_type
169 : INTEGER :: maxl = -1
170 : REAL(KIND=dp) :: beta = 0.0_dp
171 : REAL(KIND=dp) :: weight = 0.0_dp
172 : INTEGER :: max_projector = -1
173 : REAL(KIND=dp) :: beta_radius = HUGE(dp)
174 : END TYPE pao_potential_type
175 :
176 : ! **************************************************************************************************
177 : !> \brief Holds information about a PAO descriptor
178 : ! **************************************************************************************************
179 : TYPE pao_descriptor_type
180 : REAL(KIND=dp) :: beta = 0.0_dp
181 : REAL(KIND=dp) :: beta_radius = HUGE(dp)
182 : REAL(KIND=dp) :: weight = 0.0_dp
183 : REAL(KIND=dp) :: screening = 0.0_dp
184 : REAL(KIND=dp) :: screening_radius = HUGE(dp)
185 : END TYPE pao_descriptor_type
186 :
187 : ! **************************************************************************************************
188 : !> \brief Provides all information about a quickstep kind
189 : ! **************************************************************************************************
190 : TYPE qs_kind_type
191 : CHARACTER(LEN=default_string_length) :: name = ""
192 : CHARACTER(LEN=2) :: element_symbol = ""
193 : INTEGER :: natom = -1
194 : TYPE(all_potential_type), POINTER :: all_potential => Null()
195 : TYPE(local_potential_type), POINTER :: tnadd_potential => Null()
196 : TYPE(gth_potential_type), POINTER :: gth_potential => Null()
197 : TYPE(sgp_potential_type), POINTER :: sgp_potential => Null()
198 : TYPE(semi_empirical_type), POINTER :: se_parameter => Null()
199 : TYPE(qs_dftb_atom_type), POINTER :: dftb_parameter => Null()
200 : TYPE(xtb_atom_type), POINTER :: xtb_parameter => Null()
201 : !
202 : TYPE(atom_upfpot_type), POINTER :: upf_potential => Null()
203 : TYPE(cneo_potential_type), POINTER :: cneo_potential => Null()
204 : !
205 : TYPE(basis_set_container_type), &
206 : DIMENSION(20) :: basis_sets = basis_set_container_type()
207 : ! Atomic radii
208 : REAL(KIND=dp) :: covalent_radius = 0.0_dp
209 : REAL(KIND=dp) :: vdw_radius = 0.0_dp
210 : ! GAPW specific data
211 : TYPE(paw_proj_set_type), POINTER :: paw_proj_set => Null()
212 : REAL(KIND=dp) :: hard_radius = 0.8_dp*bohr ! for hard and soft exp
213 : REAL(KIND=dp) :: hard0_radius = 0.8_dp*bohr ! for hard exp of rho0
214 : REAL(KIND=dp) :: max_rad_local = 13.2_dp*bohr ! max GTO radius used in GAPW
215 : LOGICAL :: paw_atom = .FALSE. ! needs atomic rho1
216 : LOGICAL :: gpw_type_forced = .FALSE. ! gpw atom even if with hard exponents
217 : !
218 : LOGICAL :: ghost = .FALSE.
219 : LOGICAL :: monovalent = .FALSE.
220 : LOGICAL :: floating = .FALSE.
221 : INTEGER :: lmax_dftb = -1
222 : REAL(KIND=dp) :: dudq_dftb3 = 0.0_dp
223 : REAL(KIND=dp) :: magnetization = 0.0_dp
224 : INTEGER, DIMENSION(:, :), POINTER :: addel => Null()
225 : INTEGER, DIMENSION(:, :), POINTER :: laddel => Null()
226 : INTEGER, DIMENSION(:, :), POINTER :: naddel => Null()
227 : TYPE(harmonics_atom_type), POINTER :: harmonics => Null()
228 : TYPE(grid_atom_type), POINTER :: grid_atom => Null()
229 : INTEGER :: ngrid_rad = 50
230 : INTEGER :: ngrid_ang = 50
231 : INTEGER :: lmax_rho0 = 0
232 : INTEGER :: mao = -1
233 : INTEGER, DIMENSION(:), POINTER :: elec_conf => Null() ! used to set up the initial atomic guess
234 : LOGICAL :: bs_occupation = .FALSE.
235 : TYPE(dft_plus_u_type), POINTER :: dft_plus_u => Null()
236 : LOGICAL :: no_optimize = .TRUE.
237 : !
238 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: nlcc_pot => Null()
239 : !
240 : TYPE(qs_atom_dispersion_type), POINTER :: dispersion => Null()
241 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: reltmat => Null()
242 : INTEGER :: pao_basis_size = -1
243 : CHARACTER(LEN=default_path_length) :: pao_model_file = ""
244 : TYPE(pao_potential_type), DIMENSION(:), POINTER :: pao_potentials => Null()
245 : TYPE(pao_descriptor_type), DIMENSION(:), POINTER :: pao_descriptors => Null()
246 : END TYPE qs_kind_type
247 :
248 : ! Public subroutines
249 :
250 : PUBLIC :: check_qs_kind_set, &
251 : deallocate_qs_kind_set, &
252 : get_qs_kind, &
253 : get_qs_kind_set, &
254 : has_nlcc, &
255 : init_qs_kind_set, &
256 : init_gapw_basis_set, &
257 : init_gapw_nlcc, &
258 : create_qs_kind_set, &
259 : set_qs_kind, &
260 : write_qs_kind_set, &
261 : write_gto_basis_sets, &
262 : init_atom_electronic_state, set_pseudo_state, &
263 : init_cneo_basis_set
264 :
265 : ! Public data types
266 : PUBLIC :: qs_kind_type, pao_potential_type, pao_descriptor_type
267 :
268 : CONTAINS
269 :
270 : ! **************************************************************************************************
271 : !> \brief Destructor routine for a set of qs kinds
272 : !> \param qs_kind_set ...
273 : !> \date 02.01.2002
274 : !> \author Matthias Krack (MK)
275 : !> \version 2.0
276 : ! **************************************************************************************************
277 9528 : SUBROUTINE deallocate_qs_kind_set(qs_kind_set)
278 :
279 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
280 :
281 : INTEGER :: ikind, nkind
282 :
283 9528 : IF (ASSOCIATED(qs_kind_set)) THEN
284 :
285 9528 : nkind = SIZE(qs_kind_set)
286 :
287 27131 : DO ikind = 1, nkind
288 17603 : IF (ASSOCIATED(qs_kind_set(ikind)%all_potential)) THEN
289 7024 : CALL deallocate_potential(qs_kind_set(ikind)%all_potential)
290 : END IF
291 17603 : IF (ASSOCIATED(qs_kind_set(ikind)%tnadd_potential)) THEN
292 24 : CALL deallocate_potential(qs_kind_set(ikind)%tnadd_potential)
293 : END IF
294 17603 : IF (ASSOCIATED(qs_kind_set(ikind)%gth_potential)) THEN
295 10289 : CALL deallocate_potential(qs_kind_set(ikind)%gth_potential)
296 : END IF
297 17603 : IF (ASSOCIATED(qs_kind_set(ikind)%sgp_potential)) THEN
298 106 : CALL deallocate_potential(qs_kind_set(ikind)%sgp_potential)
299 : END IF
300 17603 : IF (ASSOCIATED(qs_kind_set(ikind)%upf_potential)) THEN
301 22 : CALL atom_release_upf(qs_kind_set(ikind)%upf_potential)
302 22 : DEALLOCATE (qs_kind_set(ikind)%upf_potential)
303 : END IF
304 17603 : IF (ASSOCIATED(qs_kind_set(ikind)%cneo_potential)) THEN
305 8 : CALL deallocate_cneo_potential(qs_kind_set(ikind)%cneo_potential)
306 : END IF
307 17603 : IF (ASSOCIATED(qs_kind_set(ikind)%se_parameter)) THEN
308 2236 : CALL semi_empirical_release(qs_kind_set(ikind)%se_parameter)
309 : END IF
310 17603 : IF (ASSOCIATED(qs_kind_set(ikind)%dftb_parameter)) THEN
311 630 : CALL deallocate_dftb_atom_param(qs_kind_set(ikind)%dftb_parameter)
312 : END IF
313 17603 : IF (ASSOCIATED(qs_kind_set(ikind)%xtb_parameter)) THEN
314 2614 : CALL deallocate_xtb_atom_param(qs_kind_set(ikind)%xtb_parameter)
315 : END IF
316 17603 : IF (ASSOCIATED(qs_kind_set(ikind)%paw_proj_set)) THEN
317 2548 : CALL deallocate_paw_proj_set(qs_kind_set(ikind)%paw_proj_set)
318 : END IF
319 17603 : IF (ASSOCIATED(qs_kind_set(ikind)%harmonics)) THEN
320 3036 : CALL deallocate_harmonics_atom(qs_kind_set(ikind)%harmonics)
321 : END IF
322 17603 : IF (ASSOCIATED(qs_kind_set(ikind)%grid_atom)) THEN
323 3056 : CALL deallocate_grid_atom(qs_kind_set(ikind)%grid_atom)
324 : END IF
325 17603 : IF (ASSOCIATED(qs_kind_set(ikind)%elec_conf)) THEN
326 17159 : DEALLOCATE (qs_kind_set(ikind)%elec_conf)
327 : END IF
328 :
329 17603 : IF (ASSOCIATED(qs_kind_set(ikind)%dft_plus_u)) THEN
330 90 : IF (ASSOCIATED(qs_kind_set(ikind)%dft_plus_u%orbitals)) THEN
331 4 : DEALLOCATE (qs_kind_set(ikind)%dft_plus_u%orbitals)
332 : END IF
333 90 : IF (ASSOCIATED(qs_kind_set(ikind)%dft_plus_u%nelec)) THEN
334 4 : DEALLOCATE (qs_kind_set(ikind)%dft_plus_u%nelec)
335 : END IF
336 90 : IF (ASSOCIATED(qs_kind_set(ikind)%dft_plus_u%ao_coef)) THEN
337 8 : DEALLOCATE (qs_kind_set(ikind)%dft_plus_u%ao_coef)
338 : END IF
339 90 : DEALLOCATE (qs_kind_set(ikind)%dft_plus_u)
340 : END IF
341 :
342 17603 : IF (ASSOCIATED(qs_kind_set(ikind)%nlcc_pot)) THEN
343 28 : DEALLOCATE (qs_kind_set(ikind)%nlcc_pot)
344 : END IF
345 :
346 17603 : IF (ASSOCIATED(qs_kind_set(ikind)%dispersion)) THEN
347 2790 : DEALLOCATE (qs_kind_set(ikind)%dispersion)
348 : END IF
349 17603 : IF (ASSOCIATED(qs_kind_set(ikind)%addel)) THEN
350 66 : DEALLOCATE (qs_kind_set(ikind)%addel)
351 : END IF
352 17603 : IF (ASSOCIATED(qs_kind_set(ikind)%naddel)) THEN
353 66 : DEALLOCATE (qs_kind_set(ikind)%naddel)
354 : END IF
355 17603 : IF (ASSOCIATED(qs_kind_set(ikind)%laddel)) THEN
356 66 : DEALLOCATE (qs_kind_set(ikind)%laddel)
357 : END IF
358 17603 : IF (ASSOCIATED(qs_kind_set(ikind)%reltmat)) THEN
359 28 : DEALLOCATE (qs_kind_set(ikind)%reltmat)
360 : END IF
361 :
362 17603 : IF (ASSOCIATED(qs_kind_set(ikind)%pao_potentials)) THEN
363 12045 : DEALLOCATE (qs_kind_set(ikind)%pao_potentials)
364 : END IF
365 17603 : IF (ASSOCIATED(qs_kind_set(ikind)%pao_descriptors)) THEN
366 12045 : DEALLOCATE (qs_kind_set(ikind)%pao_descriptors)
367 : END IF
368 :
369 27131 : CALL remove_basis_set_container(qs_kind_set(ikind)%basis_sets)
370 :
371 : END DO
372 9528 : DEALLOCATE (qs_kind_set)
373 : ELSE
374 : CALL cp_abort(__LOCATION__, &
375 : "The pointer qs_kind_set is not associated and "// &
376 0 : "cannot be deallocated")
377 : END IF
378 :
379 9528 : END SUBROUTINE deallocate_qs_kind_set
380 :
381 : ! **************************************************************************************************
382 : !> \brief Get attributes of an atomic kind.
383 : !> \param qs_kind ...
384 : !> \param basis_set ...
385 : !> \param basis_type ...
386 : !> \param ncgf ...
387 : !> \param nsgf ...
388 : !> \param all_potential ...
389 : !> \param tnadd_potential ...
390 : !> \param gth_potential ...
391 : !> \param sgp_potential ...
392 : !> \param upf_potential ...
393 : !> \param cneo_potential ...
394 : !> \param se_parameter ...
395 : !> \param dftb_parameter ...
396 : !> \param xtb_parameter ...
397 : !> \param dftb3_param ...
398 : !> \param zatom ...
399 : !> \param zeff ...
400 : !> \param elec_conf ...
401 : !> \param mao ...
402 : !> \param lmax_dftb ...
403 : !> \param alpha_core_charge ...
404 : !> \param ccore_charge ...
405 : !> \param core_charge ...
406 : !> \param core_charge_radius ...
407 : !> \param paw_proj_set ...
408 : !> \param paw_atom ...
409 : !> \param hard_radius ...
410 : !> \param hard0_radius ...
411 : !> \param max_rad_local ...
412 : !> \param covalent_radius ...
413 : !> \param vdw_radius ...
414 : !> \param gpw_type_forced ...
415 : !> \param harmonics ...
416 : !> \param max_iso_not0 ...
417 : !> \param max_s_harm ...
418 : !> \param grid_atom ...
419 : !> \param ngrid_ang ...
420 : !> \param ngrid_rad ...
421 : !> \param lmax_rho0 ...
422 : !> \param dft_plus_u_atom ...
423 : !> \param l_of_dft_plus_u ...
424 : !> \param n_of_dft_plus_u ...
425 : !> \param u_minus_j ...
426 : !> \param hund_j ...
427 : !> \param U_of_dft_plus_u ...
428 : !> \param J_of_dft_plus_u ...
429 : !> \param alpha_of_dft_plus_u ...
430 : !> \param beta_of_dft_plus_u ...
431 : !> \param J0_of_dft_plus_u ...
432 : !> \param occupation_of_dft_plus_u ...
433 : !> \param dispersion ...
434 : !> \param bs_occupation ...
435 : !> \param magnetization ...
436 : !> \param no_optimize ...
437 : !> \param addel ...
438 : !> \param laddel ...
439 : !> \param naddel ...
440 : !> \param orbitals ...
441 : !> \param max_scf ...
442 : !> \param eps_scf ...
443 : !> \param smear ...
444 : !> \param u_ramping ...
445 : !> \param u_minus_j_target ...
446 : !> \param eps_u_ramping ...
447 : !> \param proj_shell_charge ...
448 : !> \param lr_atom ...
449 : !> \param do_mtlr ...
450 : !> \param u_j_loop ...
451 : !> \param ao_coef ...
452 : !> \param init_u_ramping_each_scf ...
453 : !> \param reltmat ...
454 : !> \param ghost ...
455 : !> \param monovalent ...
456 : !> \param floating ...
457 : !> \param name ...
458 : !> \param element_symbol ...
459 : !> \param pao_basis_size ...
460 : !> \param pao_model_file ...
461 : !> \param pao_potentials ...
462 : !> \param pao_descriptors ...
463 : !> \param nelec ...
464 : ! **************************************************************************************************
465 74692827 : SUBROUTINE get_qs_kind(qs_kind, &
466 : basis_set, basis_type, ncgf, nsgf, &
467 : all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, &
468 : cneo_potential, se_parameter, dftb_parameter, xtb_parameter, &
469 : dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, &
470 : alpha_core_charge, ccore_charge, core_charge, core_charge_radius, &
471 : paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, &
472 : covalent_radius, vdw_radius, &
473 : gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, &
474 : ngrid_ang, ngrid_rad, lmax_rho0, &
475 : dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, &
476 : u_minus_j, hund_j, U_of_dft_plus_u, J_of_dft_plus_u, &
477 : alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, &
478 : bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, &
479 : max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, proj_shell_charge, &
480 : lr_atom, do_mtlr, u_j_loop, ao_coef, &
481 : init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, &
482 : pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
483 :
484 : TYPE(qs_kind_type) :: qs_kind
485 : TYPE(gto_basis_set_type), OPTIONAL, POINTER :: basis_set
486 : CHARACTER(len=*), OPTIONAL :: basis_type
487 : INTEGER, INTENT(OUT), OPTIONAL :: ncgf, nsgf
488 : TYPE(all_potential_type), OPTIONAL, POINTER :: all_potential
489 : TYPE(local_potential_type), OPTIONAL, POINTER :: tnadd_potential
490 : TYPE(gth_potential_type), OPTIONAL, POINTER :: gth_potential
491 : TYPE(sgp_potential_type), OPTIONAL, POINTER :: sgp_potential
492 : TYPE(atom_upfpot_type), OPTIONAL, POINTER :: upf_potential
493 : TYPE(cneo_potential_type), OPTIONAL, POINTER :: cneo_potential
494 : TYPE(semi_empirical_type), OPTIONAL, POINTER :: se_parameter
495 : TYPE(qs_dftb_atom_type), OPTIONAL, POINTER :: dftb_parameter
496 : TYPE(xtb_atom_type), OPTIONAL, POINTER :: xtb_parameter
497 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: dftb3_param
498 : INTEGER, INTENT(OUT), OPTIONAL :: zatom
499 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: zeff
500 : INTEGER, DIMENSION(:), OPTIONAL, POINTER :: elec_conf
501 : INTEGER, INTENT(OUT), OPTIONAL :: mao, lmax_dftb
502 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: alpha_core_charge, ccore_charge, &
503 : core_charge, core_charge_radius
504 : TYPE(paw_proj_set_type), OPTIONAL, POINTER :: paw_proj_set
505 : LOGICAL, INTENT(OUT), OPTIONAL :: paw_atom
506 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: hard_radius, hard0_radius, &
507 : max_rad_local, covalent_radius, &
508 : vdw_radius
509 : LOGICAL, INTENT(OUT), OPTIONAL :: gpw_type_forced
510 : TYPE(harmonics_atom_type), OPTIONAL, POINTER :: harmonics
511 : INTEGER, INTENT(OUT), OPTIONAL :: max_iso_not0, max_s_harm
512 : TYPE(grid_atom_type), OPTIONAL, POINTER :: grid_atom
513 : INTEGER, INTENT(OUT), OPTIONAL :: ngrid_ang, ngrid_rad, lmax_rho0
514 : LOGICAL, INTENT(OUT), OPTIONAL :: dft_plus_u_atom
515 : INTEGER, INTENT(OUT), OPTIONAL :: l_of_dft_plus_u, n_of_dft_plus_u
516 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: u_minus_j, hund_j, U_of_dft_plus_u, J_of_dft_plus_u, &
517 : alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u
518 : TYPE(qs_atom_dispersion_type), OPTIONAL, POINTER :: dispersion
519 : LOGICAL, INTENT(OUT), OPTIONAL :: bs_occupation
520 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: magnetization
521 : LOGICAL, INTENT(OUT), OPTIONAL :: no_optimize
522 : INTEGER, DIMENSION(:, :), OPTIONAL, POINTER :: addel, laddel, naddel
523 : INTEGER, DIMENSION(:), OPTIONAL, POINTER :: orbitals
524 : INTEGER, OPTIONAL :: max_scf
525 : REAL(KIND=dp), OPTIONAL :: eps_scf
526 : LOGICAL, OPTIONAL :: smear
527 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: u_ramping, u_minus_j_target, &
528 : eps_u_ramping, proj_shell_charge(4)
529 : INTEGER, OPTIONAL :: lr_atom
530 : LOGICAL, INTENT(OUT), OPTIONAL :: do_mtlr, u_j_loop
531 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: ao_coef
532 : LOGICAL, OPTIONAL :: init_u_ramping_each_scf
533 : REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: reltmat
534 : LOGICAL, OPTIONAL :: ghost
535 : LOGICAL, INTENT(OUT), OPTIONAL :: monovalent
536 : LOGICAL, OPTIONAL :: floating
537 : CHARACTER(LEN=default_string_length), &
538 : INTENT(OUT), OPTIONAL :: name
539 : CHARACTER(LEN=2), INTENT(OUT), OPTIONAL :: element_symbol
540 : INTEGER, INTENT(OUT), OPTIONAL :: pao_basis_size
541 : CHARACTER(LEN=default_path_length), INTENT(OUT), &
542 : OPTIONAL :: pao_model_file
543 : TYPE(pao_potential_type), DIMENSION(:), OPTIONAL, &
544 : POINTER :: pao_potentials
545 : TYPE(pao_descriptor_type), DIMENSION(:), &
546 : OPTIONAL, POINTER :: pao_descriptors
547 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: nelec
548 :
549 : CHARACTER(LEN=default_string_length) :: my_basis_type
550 : INTEGER :: l
551 : LOGICAL :: found
552 : TYPE(gto_basis_set_type), POINTER :: tmp_basis_set
553 :
554 : ! Retrieve basis set from the kind container
555 74692827 : IF (PRESENT(basis_type)) THEN
556 19350415 : my_basis_type = basis_type
557 : ELSE
558 55342412 : my_basis_type = "ORB"
559 : END IF
560 :
561 74692827 : IF (PRESENT(basis_set)) THEN
562 : CALL get_basis_from_container(qs_kind%basis_sets, basis_set=basis_set, &
563 20395240 : basis_type=my_basis_type)
564 : END IF
565 :
566 74692827 : IF (PRESENT(ncgf)) THEN
567 : CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis_set, &
568 960 : basis_type=my_basis_type)
569 960 : IF (ASSOCIATED(tmp_basis_set)) THEN
570 960 : CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, ncgf=ncgf)
571 0 : ELSE IF (ASSOCIATED(qs_kind%dftb_parameter)) THEN
572 0 : l = qs_kind%dftb_parameter%lmax
573 0 : ncgf = ((l + 1)*(l + 2)*(l + 3))/6
574 : ELSE
575 0 : ncgf = 0
576 : END IF
577 : END IF
578 :
579 74692827 : IF (PRESENT(nsgf)) THEN
580 : CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis_set, &
581 509129 : basis_type=my_basis_type)
582 509129 : IF (ASSOCIATED(tmp_basis_set)) THEN
583 395053 : CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, nsgf=nsgf)
584 114076 : ELSE IF (ASSOCIATED(qs_kind%dftb_parameter)) THEN
585 114072 : nsgf = qs_kind%dftb_parameter%natorb
586 : ELSE
587 4 : nsgf = 0
588 : END IF
589 : END IF
590 :
591 74692827 : IF (PRESENT(all_potential)) all_potential => qs_kind%all_potential
592 74692827 : IF (PRESENT(tnadd_potential)) tnadd_potential => qs_kind%tnadd_potential
593 74692827 : IF (PRESENT(gth_potential)) gth_potential => qs_kind%gth_potential
594 74692827 : IF (PRESENT(sgp_potential)) sgp_potential => qs_kind%sgp_potential
595 74692827 : IF (PRESENT(upf_potential)) upf_potential => qs_kind%upf_potential
596 74692827 : IF (PRESENT(cneo_potential)) cneo_potential => qs_kind%cneo_potential
597 74692827 : IF (PRESENT(se_parameter)) se_parameter => qs_kind%se_parameter
598 74692827 : IF (PRESENT(dftb_parameter)) dftb_parameter => qs_kind%dftb_parameter
599 74692827 : IF (PRESENT(xtb_parameter)) xtb_parameter => qs_kind%xtb_parameter
600 74692827 : IF (PRESENT(element_symbol)) element_symbol = qs_kind%element_symbol
601 74692827 : IF (PRESENT(name)) name = qs_kind%name
602 74692827 : IF (PRESENT(dftb3_param)) dftb3_param = qs_kind%dudq_dftb3
603 74692827 : IF (PRESENT(elec_conf)) elec_conf => qs_kind%elec_conf
604 74692827 : IF (PRESENT(alpha_core_charge)) THEN
605 253955 : IF (ASSOCIATED(qs_kind%all_potential)) THEN
606 : CALL get_potential(potential=qs_kind%all_potential, &
607 67376 : alpha_core_charge=alpha_core_charge)
608 186579 : ELSE IF (ASSOCIATED(qs_kind%gth_potential)) THEN
609 : CALL get_potential(potential=qs_kind%gth_potential, &
610 183801 : alpha_core_charge=alpha_core_charge)
611 2778 : ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
612 : CALL get_potential(potential=qs_kind%sgp_potential, &
613 1346 : alpha_core_charge=alpha_core_charge)
614 1432 : ELSE IF (ASSOCIATED(qs_kind%cneo_potential)) THEN
615 0 : CPABORT("CNEO ALPHA CORE CHARGE NOT AVAILABLE")
616 : ELSE
617 1432 : alpha_core_charge = 1.0_dp
618 : END IF
619 : END IF
620 74692827 : IF (PRESENT(ccore_charge)) THEN
621 102759 : IF (ASSOCIATED(qs_kind%all_potential)) THEN
622 : CALL get_potential(potential=qs_kind%all_potential, &
623 16578 : ccore_charge=ccore_charge)
624 86181 : ELSE IF (ASSOCIATED(qs_kind%gth_potential)) THEN
625 : CALL get_potential(potential=qs_kind%gth_potential, &
626 84703 : ccore_charge=ccore_charge)
627 1478 : ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
628 : CALL get_potential(potential=qs_kind%sgp_potential, &
629 654 : ccore_charge=ccore_charge)
630 824 : ELSE IF (ASSOCIATED(qs_kind%upf_potential)) THEN
631 0 : CPABORT("UPF CCORE CHARGE NOT AVAILABLE")
632 824 : ELSE IF (ASSOCIATED(qs_kind%cneo_potential)) THEN
633 0 : CPABORT("CNEO CCORE CHARGE NOT AVAILABLE")
634 : ELSE
635 824 : ccore_charge = 0.0_dp
636 : END IF
637 : END IF
638 74692827 : IF (PRESENT(core_charge_radius)) THEN
639 104739 : IF (ASSOCIATED(qs_kind%all_potential)) THEN
640 : CALL get_potential(potential=qs_kind%all_potential, &
641 45850 : core_charge_radius=core_charge_radius)
642 58889 : ELSE IF (ASSOCIATED(qs_kind%gth_potential)) THEN
643 : CALL get_potential(potential=qs_kind%gth_potential, &
644 58099 : core_charge_radius=core_charge_radius)
645 790 : ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
646 : CALL get_potential(potential=qs_kind%sgp_potential, &
647 396 : core_charge_radius=core_charge_radius)
648 394 : ELSE IF (ASSOCIATED(qs_kind%upf_potential)) THEN
649 0 : CPABORT("UPF CORE CHARGE RADIUS NOT AVAILABLE")
650 394 : ELSE IF (ASSOCIATED(qs_kind%cneo_potential)) THEN
651 0 : CPABORT("CNEO CORE CHARGE RADIUS NOT AVAILABLE")
652 : ELSE
653 394 : core_charge_radius = 0.0_dp
654 : END IF
655 : END IF
656 74692827 : IF (PRESENT(core_charge)) THEN
657 51407 : IF (ASSOCIATED(qs_kind%all_potential)) THEN
658 : CALL get_potential(potential=qs_kind%all_potential, &
659 1138 : zeff=core_charge)
660 50269 : ELSE IF (ASSOCIATED(qs_kind%gth_potential)) THEN
661 : CALL get_potential(potential=qs_kind%gth_potential, &
662 50269 : zeff=core_charge)
663 0 : ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
664 : CALL get_potential(potential=qs_kind%sgp_potential, &
665 0 : zeff=core_charge)
666 0 : ELSE IF (ASSOCIATED(qs_kind%upf_potential)) THEN
667 0 : CPABORT("UPF CORE CHARGE NOT AVAILABLE")
668 0 : ELSE IF (ASSOCIATED(qs_kind%cneo_potential)) THEN
669 : CALL get_cneo_potential(potential=qs_kind%cneo_potential, &
670 0 : zeff=core_charge)
671 : ELSE
672 0 : core_charge = 0.0_dp
673 : END IF
674 : END IF
675 :
676 74692827 : IF (PRESENT(zatom)) THEN
677 : ! Retrieve information on element
678 310566 : CALL get_ptable_info(qs_kind%element_symbol, ielement=zatom, found=found)
679 310566 : CPASSERT(found)
680 : END IF
681 :
682 74692827 : IF (PRESENT(zeff)) THEN
683 350501 : IF (ASSOCIATED(qs_kind%all_potential)) THEN
684 91967 : CALL get_potential(potential=qs_kind%all_potential, zeff=zeff)
685 258534 : ELSE IF (ASSOCIATED(qs_kind%gth_potential)) THEN
686 254718 : CALL get_potential(potential=qs_kind%gth_potential, zeff=zeff)
687 3816 : ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
688 2424 : CALL get_potential(potential=qs_kind%sgp_potential, zeff=zeff)
689 1392 : ELSE IF (ASSOCIATED(qs_kind%upf_potential)) THEN
690 73 : zeff = qs_kind%upf_potential%zion
691 1319 : ELSE IF (ASSOCIATED(qs_kind%cneo_potential)) THEN
692 143 : CALL get_cneo_potential(potential=qs_kind%cneo_potential, zeff=zeff)
693 : ELSE
694 1176 : zeff = 0.0_dp
695 : END IF
696 : END IF
697 :
698 74692827 : IF (PRESENT(covalent_radius)) covalent_radius = qs_kind%covalent_radius
699 74692827 : IF (PRESENT(vdw_radius)) vdw_radius = qs_kind%vdw_radius
700 :
701 74692827 : IF (PRESENT(paw_proj_set)) paw_proj_set => qs_kind%paw_proj_set
702 74692827 : IF (PRESENT(paw_atom)) paw_atom = qs_kind%paw_atom
703 74692827 : IF (PRESENT(gpw_type_forced)) gpw_type_forced = qs_kind%gpw_type_forced
704 74692827 : IF (PRESENT(hard_radius)) hard_radius = qs_kind%hard_radius
705 74692827 : IF (PRESENT(hard0_radius)) hard0_radius = qs_kind%hard0_radius
706 74692827 : IF (PRESENT(max_rad_local)) max_rad_local = qs_kind%max_rad_local
707 74692827 : IF (PRESENT(harmonics)) harmonics => qs_kind%harmonics
708 74692827 : IF (PRESENT(max_s_harm)) THEN
709 9379097 : IF (ASSOCIATED(qs_kind%harmonics)) THEN
710 488587 : max_s_harm = qs_kind%harmonics%max_s_harm
711 : ELSE
712 8890510 : max_s_harm = 0
713 : END IF
714 : END IF
715 74692827 : IF (PRESENT(max_iso_not0)) THEN
716 9434455 : IF (ASSOCIATED(qs_kind%harmonics)) THEN
717 543945 : max_iso_not0 = qs_kind%harmonics%max_iso_not0
718 : ELSE
719 8890510 : max_iso_not0 = 0
720 : END IF
721 : END IF
722 74692827 : IF (PRESENT(grid_atom)) grid_atom => qs_kind%grid_atom
723 74692827 : IF (PRESENT(ngrid_ang)) ngrid_ang = qs_kind%ngrid_ang
724 74692827 : IF (PRESENT(ngrid_rad)) ngrid_rad = qs_kind%ngrid_rad
725 74692827 : IF (PRESENT(lmax_rho0)) lmax_rho0 = qs_kind%lmax_rho0
726 74692827 : IF (PRESENT(ghost)) ghost = qs_kind%ghost
727 74692827 : IF (PRESENT(monovalent)) monovalent = qs_kind%monovalent
728 74692827 : IF (PRESENT(floating)) floating = qs_kind%floating
729 74692827 : IF (PRESENT(dft_plus_u_atom)) dft_plus_u_atom = ASSOCIATED(qs_kind%dft_plus_u)
730 74692827 : IF (PRESENT(do_mtlr)) THEN
731 9379097 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
732 29402 : do_mtlr = qs_kind%dft_plus_u%do_mtlr
733 : ELSE
734 9349695 : do_mtlr = .FALSE.
735 : END IF
736 : END IF
737 74692827 : IF (PRESENT(l_of_dft_plus_u)) THEN
738 209398 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
739 106412 : l_of_dft_plus_u = qs_kind%dft_plus_u%l
740 : ELSE
741 102986 : l_of_dft_plus_u = -1
742 : END IF
743 : END IF
744 74692827 : IF (PRESENT(n_of_dft_plus_u)) THEN
745 30 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
746 0 : n_of_dft_plus_u = qs_kind%dft_plus_u%n
747 : ELSE
748 30 : n_of_dft_plus_u = -1
749 : END IF
750 : END IF
751 74692827 : IF (PRESENT(u_minus_j)) THEN
752 210144 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
753 106800 : u_minus_j = qs_kind%dft_plus_u%u_minus_j
754 : ELSE
755 103344 : u_minus_j = 0.0_dp
756 : END IF
757 : END IF
758 74692827 : IF (PRESENT(u_minus_j_target)) THEN
759 210378 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
760 106902 : u_minus_j_target = qs_kind%dft_plus_u%u_minus_j_target
761 : ELSE
762 103476 : u_minus_j_target = 0.0_dp
763 : END IF
764 : END IF
765 74692827 : IF (PRESENT(hund_j)) THEN
766 13104 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
767 8878 : hund_j = qs_kind%dft_plus_u%hund_j
768 : ELSE
769 4226 : hund_j = 0.0_dp
770 : END IF
771 : END IF
772 74692827 : IF (PRESENT(U_of_dft_plus_u)) THEN
773 30 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
774 0 : U_of_dft_plus_u = qs_kind%dft_plus_u%U
775 : ELSE
776 30 : U_of_dft_plus_u = 0.0_dp
777 : END IF
778 : END IF
779 74692827 : IF (PRESENT(J_of_dft_plus_u)) THEN
780 30 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
781 0 : J_of_dft_plus_u = qs_kind%dft_plus_u%J
782 : ELSE
783 30 : J_of_dft_plus_u = 0.0_dp
784 : END IF
785 : END IF
786 74692827 : IF (PRESENT(alpha_of_dft_plus_u)) THEN
787 30 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
788 0 : alpha_of_dft_plus_u = qs_kind%dft_plus_u%alpha
789 : ELSE
790 30 : alpha_of_dft_plus_u = 0.0_dp
791 : END IF
792 : END IF
793 74692827 : IF (PRESENT(beta_of_dft_plus_u)) THEN
794 30 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
795 0 : beta_of_dft_plus_u = qs_kind%dft_plus_u%beta
796 : ELSE
797 30 : beta_of_dft_plus_u = 0.0_dp
798 : END IF
799 : END IF
800 74692827 : IF (PRESENT(J0_of_dft_plus_u)) THEN
801 30 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
802 0 : J0_of_dft_plus_u = qs_kind%dft_plus_u%J0
803 : ELSE
804 30 : J0_of_dft_plus_u = 0.0_dp
805 : END IF
806 : END IF
807 74692827 : IF (PRESENT(occupation_of_dft_plus_u)) THEN
808 30 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
809 0 : occupation_of_dft_plus_u = qs_kind%dft_plus_u%occupation
810 : ELSE
811 30 : occupation_of_dft_plus_u = -1.0_dp
812 : END IF
813 : END IF
814 :
815 74692827 : IF (PRESENT(init_u_ramping_each_scf)) THEN
816 1260 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
817 730 : init_u_ramping_each_scf = qs_kind%dft_plus_u%init_u_ramping_each_scf
818 : ELSE
819 530 : init_u_ramping_each_scf = .FALSE.
820 : END IF
821 : END IF
822 74692827 : IF (PRESENT(u_ramping)) THEN
823 211384 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
824 107520 : u_ramping = qs_kind%dft_plus_u%u_ramping
825 : ELSE
826 103864 : u_ramping = 0.0_dp
827 : END IF
828 : END IF
829 74692827 : IF (PRESENT(eps_u_ramping)) THEN
830 210108 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
831 106782 : eps_u_ramping = qs_kind%dft_plus_u%eps_u_ramping
832 : ELSE
833 103326 : eps_u_ramping = 1.0E-5_dp
834 : END IF
835 : END IF
836 74692827 : IF (PRESENT(lr_atom)) THEN
837 0 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
838 0 : lr_atom = qs_kind%dft_plus_u%lr_atom
839 : ELSE
840 0 : lr_atom = -1
841 : END IF
842 : END IF
843 74692827 : IF (PRESENT(u_j_loop)) THEN
844 0 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
845 0 : u_j_loop = qs_kind%dft_plus_u%u_j_loop
846 : ELSE
847 0 : u_j_loop = .FALSE.
848 : END IF
849 : END IF
850 74692827 : IF (PRESENT(proj_shell_charge)) THEN
851 0 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
852 0 : proj_shell_charge = qs_kind%dft_plus_u%proj_shell_charge
853 : ELSE
854 0 : proj_shell_charge = 0.0_dp
855 : END IF
856 : END IF
857 74692827 : IF (PRESENT(nelec)) THEN
858 197464 : NULLIFY (nelec)
859 197464 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
860 98732 : IF (ASSOCIATED(qs_kind%dft_plus_u%nelec)) THEN
861 0 : nelec => qs_kind%dft_plus_u%nelec
862 : END IF
863 : END IF
864 : END IF
865 74692827 : IF (PRESENT(orbitals)) THEN
866 210532 : NULLIFY (orbitals)
867 210532 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
868 107592 : IF (ASSOCIATED(qs_kind%dft_plus_u%orbitals)) THEN
869 204 : orbitals => qs_kind%dft_plus_u%orbitals
870 : END IF
871 : END IF
872 : END IF
873 74692827 : IF (PRESENT(eps_scf)) THEN
874 206360 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
875 104908 : eps_scf = qs_kind%dft_plus_u%eps_scf
876 : ELSE
877 101452 : eps_scf = 1.0E30_dp
878 : END IF
879 : END IF
880 74692827 : IF (PRESENT(max_scf)) THEN
881 206360 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
882 104908 : max_scf = qs_kind%dft_plus_u%max_scf
883 : ELSE
884 101452 : max_scf = -1
885 : END IF
886 : END IF
887 74692827 : IF (PRESENT(smear)) THEN
888 206360 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
889 104908 : smear = qs_kind%dft_plus_u%smear
890 : ELSE
891 101452 : smear = .FALSE.
892 : END IF
893 : END IF
894 74692827 : IF (PRESENT(dispersion)) dispersion => qs_kind%dispersion
895 74692827 : IF (PRESENT(bs_occupation)) bs_occupation = qs_kind%bs_occupation
896 74692827 : IF (PRESENT(addel)) addel => qs_kind%addel
897 74692827 : IF (PRESENT(laddel)) laddel => qs_kind%laddel
898 74692827 : IF (PRESENT(naddel)) naddel => qs_kind%naddel
899 :
900 74692827 : IF (PRESENT(magnetization)) magnetization = qs_kind%magnetization
901 :
902 74692827 : IF (PRESENT(no_optimize)) no_optimize = qs_kind%no_optimize
903 :
904 74692827 : IF (PRESENT(reltmat)) reltmat => qs_kind%reltmat
905 :
906 74692827 : IF (PRESENT(ao_coef)) THEN
907 2024 : NULLIFY (ao_coef)
908 2024 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
909 1544 : IF (ASSOCIATED(qs_kind%dft_plus_u%ao_coef)) THEN
910 1544 : ao_coef => qs_kind%dft_plus_u%ao_coef
911 : END IF
912 : END IF
913 : END IF
914 :
915 74692827 : IF (PRESENT(mao)) mao = qs_kind%mao
916 :
917 74692827 : IF (PRESENT(lmax_dftb)) lmax_dftb = qs_kind%lmax_dftb
918 :
919 74692827 : IF (PRESENT(pao_basis_size)) pao_basis_size = qs_kind%pao_basis_size
920 74692827 : IF (PRESENT(pao_model_file)) pao_model_file = qs_kind%pao_model_file
921 74692827 : IF (PRESENT(pao_potentials)) pao_potentials => qs_kind%pao_potentials
922 74692827 : IF (PRESENT(pao_descriptors)) pao_descriptors => qs_kind%pao_descriptors
923 74692827 : END SUBROUTINE get_qs_kind
924 :
925 : ! **************************************************************************************************
926 : !> \brief Get attributes of an atomic kind set.
927 : !> \param qs_kind_set ...
928 : !> \param all_potential_present ...
929 : !> \param tnadd_potential_present ...
930 : !> \param gth_potential_present ...
931 : !> \param sgp_potential_present ...
932 : !> \param paw_atom_present ...
933 : !> \param dft_plus_u_atom_present ...
934 : !> \param maxcgf ...
935 : !> \param maxsgf ...
936 : !> \param maxco ...
937 : !> \param maxco_proj ...
938 : !> \param maxgtops ...
939 : !> \param maxlgto ...
940 : !> \param maxlprj ...
941 : !> \param maxnset ...
942 : !> \param maxsgf_set ...
943 : !> \param ncgf ...
944 : !> \param npgf ...
945 : !> \param nset ...
946 : !> \param nsgf ...
947 : !> \param nshell ...
948 : !> \param maxpol ...
949 : !> \param maxlppl ...
950 : !> \param maxlppnl ...
951 : !> \param maxppnl ...
952 : !> \param nelectron ...
953 : !> \param maxder ...
954 : !> \param max_ngrid_rad ...
955 : !> \param max_sph_harm ...
956 : !> \param maxg_iso_not0 ...
957 : !> \param lmax_rho0 ...
958 : !> \param basis_rcut ...
959 : !> \param do_mtlr_present ...
960 : !> \param basis_type ...
961 : !> \param total_zeff_corr ... [SGh]
962 : !> \param npgf_seg total number of primitive GTOs in "segmented contraction format"
963 : !> \param cneo_potential_present ...
964 : !> \param nkind_q ...
965 : !> \param natom_q ...
966 : ! **************************************************************************************************
967 4671986 : SUBROUTINE get_qs_kind_set(qs_kind_set, &
968 : all_potential_present, tnadd_potential_present, gth_potential_present, &
969 : sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, &
970 : maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, &
971 : ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, &
972 : nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, &
973 : basis_rcut, do_mtlr_present, &
974 : basis_type, total_zeff_corr, npgf_seg, &
975 : cneo_potential_present, nkind_q, natom_q)
976 :
977 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
978 : LOGICAL, INTENT(OUT), OPTIONAL :: all_potential_present, tnadd_potential_present, &
979 : gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present
980 : INTEGER, INTENT(OUT), OPTIONAL :: maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, &
981 : maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, &
982 : maxppnl, nelectron
983 : INTEGER, INTENT(IN), OPTIONAL :: maxder
984 : INTEGER, INTENT(OUT), OPTIONAL :: max_ngrid_rad, max_sph_harm, &
985 : maxg_iso_not0, lmax_rho0
986 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: basis_rcut
987 : LOGICAL, INTENT(OUT), OPTIONAL :: do_mtlr_present
988 : CHARACTER(len=*), OPTIONAL :: basis_type
989 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: total_zeff_corr
990 : INTEGER, INTENT(OUT), OPTIONAL :: npgf_seg
991 : LOGICAL, INTENT(OUT), OPTIONAL :: cneo_potential_present
992 : INTEGER, INTENT(OUT), OPTIONAL :: nkind_q, natom_q
993 :
994 : CHARACTER(len=default_string_length) :: my_basis_type
995 : INTEGER :: ikind, imax, lmax_rho0_kind, &
996 : max_iso_not0, max_s_harm, n, &
997 : ngrid_rad, nkind, nrloc(10), &
998 : nrpot(1:15, 0:10)
999 : LOGICAL :: dft_plus_u_atom, do_mtlr, &
1000 : ecp_semi_local, paw_atom
1001 : REAL(KIND=dp) :: brcut, zeff, zeff_correction
1002 : TYPE(all_potential_type), POINTER :: all_potential
1003 : TYPE(cneo_potential_type), POINTER :: cneo_potential
1004 : TYPE(gth_potential_type), POINTER :: gth_potential
1005 : TYPE(gto_basis_set_type), POINTER :: tmp_basis_set
1006 : TYPE(local_potential_type), POINTER :: tnadd_potential
1007 : TYPE(paw_proj_set_type), POINTER :: paw_proj_set
1008 : TYPE(qs_dftb_atom_type), POINTER :: dftb_parameter
1009 : TYPE(qs_kind_type), POINTER :: qs_kind
1010 : TYPE(sgp_potential_type), POINTER :: sgp_potential
1011 :
1012 4671986 : IF (PRESENT(basis_type)) THEN
1013 4207001 : my_basis_type = basis_type
1014 : ELSE
1015 464985 : my_basis_type = "ORB"
1016 : END IF
1017 :
1018 4671986 : IF (ASSOCIATED(qs_kind_set)) THEN
1019 :
1020 4671986 : IF (PRESENT(maxcgf)) maxcgf = 0
1021 4671986 : IF (PRESENT(maxco)) maxco = 0
1022 4671986 : IF (PRESENT(maxco_proj)) maxco_proj = 0
1023 4671986 : IF (PRESENT(maxg_iso_not0)) maxg_iso_not0 = 0
1024 4671986 : IF (PRESENT(maxgtops)) maxgtops = 0
1025 4671986 : IF (PRESENT(maxlgto)) maxlgto = -1
1026 4671986 : IF (PRESENT(maxlppl)) maxlppl = -1
1027 4671986 : IF (PRESENT(maxlppnl)) maxlppnl = -1
1028 4671986 : IF (PRESENT(maxpol)) maxpol = -1
1029 4671986 : IF (PRESENT(maxlprj)) maxlprj = -1
1030 4671986 : IF (PRESENT(maxnset)) maxnset = 0
1031 4671986 : IF (PRESENT(maxppnl)) maxppnl = 0
1032 4671986 : IF (PRESENT(maxsgf)) maxsgf = 0
1033 4671986 : IF (PRESENT(maxsgf_set)) maxsgf_set = 0
1034 4671986 : IF (PRESENT(ncgf)) ncgf = 0
1035 4671986 : IF (PRESENT(nelectron)) nelectron = 0
1036 4671986 : IF (PRESENT(npgf)) npgf = 0
1037 4671986 : IF (PRESENT(nset)) nset = 0
1038 4671986 : IF (PRESENT(nsgf)) nsgf = 0
1039 4671986 : IF (PRESENT(nshell)) nshell = 0
1040 4671986 : IF (PRESENT(all_potential_present)) all_potential_present = .FALSE.
1041 4671986 : IF (PRESENT(tnadd_potential_present)) tnadd_potential_present = .FALSE.
1042 4671986 : IF (PRESENT(gth_potential_present)) gth_potential_present = .FALSE.
1043 4671986 : IF (PRESENT(sgp_potential_present)) sgp_potential_present = .FALSE.
1044 4671986 : IF (PRESENT(cneo_potential_present)) cneo_potential_present = .FALSE.
1045 4671986 : IF (PRESENT(nkind_q)) nkind_q = 0
1046 4671986 : IF (PRESENT(natom_q)) natom_q = 0
1047 4671986 : IF (PRESENT(paw_atom_present)) paw_atom_present = .FALSE.
1048 4671986 : IF (PRESENT(max_ngrid_rad)) max_ngrid_rad = 0
1049 4671986 : IF (PRESENT(max_sph_harm)) max_sph_harm = 0
1050 4671986 : IF (PRESENT(lmax_rho0)) lmax_rho0 = 0
1051 4671986 : IF (PRESENT(basis_rcut)) basis_rcut = 0.0_dp
1052 4671986 : IF (PRESENT(total_zeff_corr)) total_zeff_corr = 0.0_dp
1053 4671986 : IF (PRESENT(npgf_seg)) npgf_seg = 0
1054 4671986 : IF (PRESENT(do_mtlr_present)) do_mtlr_present = .FALSE.
1055 :
1056 4671986 : nkind = SIZE(qs_kind_set)
1057 14051083 : DO ikind = 1, nkind
1058 9379097 : qs_kind => qs_kind_set(ikind)
1059 : CALL get_qs_kind(qs_kind=qs_kind, &
1060 : all_potential=all_potential, &
1061 : tnadd_potential=tnadd_potential, &
1062 : gth_potential=gth_potential, &
1063 : sgp_potential=sgp_potential, &
1064 : cneo_potential=cneo_potential, &
1065 : paw_proj_set=paw_proj_set, &
1066 : dftb_parameter=dftb_parameter, &
1067 : ngrid_rad=ngrid_rad, &
1068 : max_s_harm=max_s_harm, &
1069 : max_iso_not0=max_iso_not0, &
1070 : paw_atom=paw_atom, &
1071 : dft_plus_u_atom=dft_plus_u_atom, &
1072 : do_mtlr=do_mtlr, &
1073 9379097 : lmax_rho0=lmax_rho0_kind)
1074 :
1075 9379097 : IF (PRESENT(maxlppl) .AND. ASSOCIATED(gth_potential)) THEN
1076 51488 : CALL get_potential(potential=gth_potential, nexp_ppl=n)
1077 51488 : maxlppl = MAX(maxlppl, 2*(n - 1))
1078 11320 : ELSE IF (PRESENT(maxlppl) .AND. ASSOCIATED(sgp_potential)) THEN
1079 393 : CALL get_potential(potential=sgp_potential, nrloc=nrloc, ecp_semi_local=ecp_semi_local)
1080 4323 : n = MAXVAL(nrloc) - 2
1081 393 : maxlppl = MAX(maxlppl, 2*(n - 1))
1082 393 : IF (ecp_semi_local) THEN
1083 321 : CALL get_potential(potential=sgp_potential, sl_lmax=imax, nrpot=nrpot)
1084 56817 : n = MAXVAL(nrpot) - 2
1085 321 : n = 2*(n - 1) + imax
1086 321 : maxlppl = MAX(maxlppl, n)
1087 : END IF
1088 : END IF
1089 :
1090 9379097 : IF (PRESENT(maxlppnl) .AND. ASSOCIATED(gth_potential)) THEN
1091 48325 : CALL get_potential(potential=gth_potential, lprj_ppnl_max=imax)
1092 48325 : maxlppnl = MAX(maxlppnl, imax)
1093 11088 : ELSE IF (PRESENT(maxlppnl) .AND. ASSOCIATED(sgp_potential)) THEN
1094 193 : CALL get_potential(potential=sgp_potential, lmax=imax)
1095 193 : maxlppnl = MAX(maxlppnl, imax)
1096 : END IF
1097 :
1098 9379097 : IF (PRESENT(maxpol) .AND. ASSOCIATED(tnadd_potential)) THEN
1099 74 : CALL get_potential(potential=tnadd_potential, npol=n)
1100 74 : maxpol = MAX(maxpol, 2*(n - 1))
1101 : END IF
1102 :
1103 9379097 : IF (PRESENT(maxco_proj) .AND. ASSOCIATED(paw_proj_set)) THEN
1104 6908 : CALL get_paw_proj_set(paw_proj_set=paw_proj_set, ncgauprj=imax)
1105 6908 : maxco_proj = MAX(maxco_proj, imax)
1106 : END IF
1107 :
1108 9379097 : IF (PRESENT(maxlprj) .AND. ASSOCIATED(paw_proj_set)) THEN
1109 6908 : CALL get_paw_proj_set(paw_proj_set=paw_proj_set, maxl=imax)
1110 6908 : maxlprj = MAX(maxlprj, imax)
1111 : END IF
1112 :
1113 9379097 : IF (PRESENT(maxppnl) .AND. ASSOCIATED(gth_potential)) THEN
1114 32946 : CALL get_potential(potential=gth_potential, nppnl=imax)
1115 32946 : maxppnl = MAX(maxppnl, imax)
1116 268 : ELSE IF (PRESENT(maxppnl) .AND. ASSOCIATED(sgp_potential)) THEN
1117 34 : CALL get_potential(potential=sgp_potential, nppnl=imax)
1118 34 : maxppnl = MAX(maxppnl, imax)
1119 : END IF
1120 :
1121 9379097 : IF (my_basis_type(1:3) == "NUC" .AND. .NOT. ASSOCIATED(cneo_potential)) THEN
1122 8700 : NULLIFY (tmp_basis_set)
1123 : ELSE
1124 : CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis_set, &
1125 9370397 : basis_type=my_basis_type)
1126 : END IF
1127 :
1128 9379097 : IF (PRESENT(maxcgf)) THEN
1129 0 : IF (ASSOCIATED(tmp_basis_set)) THEN
1130 0 : CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, ncgf=imax)
1131 0 : maxcgf = MAX(maxcgf, imax)
1132 0 : ELSE IF (ASSOCIATED(qs_kind%dftb_parameter)) THEN
1133 0 : CALL get_dftb_atom_param(dftb_parameter=dftb_parameter, lmax=imax)
1134 0 : imax = ((imax + 1)*(imax + 2)*(imax + 3))/6
1135 0 : maxcgf = MAX(maxcgf, imax)
1136 : END IF
1137 : END IF
1138 :
1139 9379097 : IF (PRESENT(maxco)) THEN
1140 8281358 : IF (ASSOCIATED(tmp_basis_set)) THEN
1141 8281250 : IF (PRESENT(maxder)) THEN
1142 : CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, &
1143 0 : maxco=imax, maxder=maxder)
1144 : ELSE
1145 8281250 : CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, maxco=imax)
1146 : END IF
1147 8281250 : maxco = MAX(maxco, imax)
1148 : END IF
1149 8281358 : IF (ASSOCIATED(gth_potential)) THEN
1150 768940 : CALL get_potential(potential=gth_potential, lprj_ppnl_max=imax)
1151 768940 : maxco = MAX(maxco, ncoset(imax))
1152 : END IF
1153 8281358 : IF (ASSOCIATED(sgp_potential)) THEN
1154 3236 : CALL get_potential(potential=sgp_potential, lmax=imax)
1155 3236 : maxco = MAX(maxco, ncoset(imax))
1156 3236 : CALL get_potential(potential=sgp_potential, sl_lmax=imax)
1157 3236 : maxco = MAX(maxco, ncoset(imax))
1158 : END IF
1159 : END IF
1160 :
1161 9379097 : IF (PRESENT(maxgtops)) THEN
1162 173108 : IF (ASSOCIATED(tmp_basis_set)) THEN
1163 173108 : CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, maxso=imax, nset=n)
1164 173108 : maxgtops = MAX(maxgtops, n*imax)
1165 : END IF
1166 : END IF
1167 :
1168 9379097 : IF (PRESENT(maxlgto)) THEN
1169 7796016 : IF (ASSOCIATED(tmp_basis_set)) THEN
1170 7757860 : CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, maxl=imax)
1171 7757860 : maxlgto = MAX(maxlgto, imax)
1172 38156 : ELSE IF (ASSOCIATED(qs_kind%dftb_parameter)) THEN
1173 5543 : CALL get_dftb_atom_param(dftb_parameter=dftb_parameter, lmax=imax)
1174 5543 : maxlgto = MAX(maxlgto, imax)
1175 : END IF
1176 : END IF
1177 :
1178 9379097 : IF (PRESENT(maxnset)) THEN
1179 82827 : IF (ASSOCIATED(tmp_basis_set)) THEN
1180 82815 : CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, nset=n)
1181 82815 : maxnset = MAX(maxnset, n)
1182 : END IF
1183 : END IF
1184 :
1185 9379097 : IF (PRESENT(maxsgf)) THEN
1186 8127066 : IF (ASSOCIATED(tmp_basis_set)) THEN
1187 8127030 : CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, nsgf=imax)
1188 8127030 : maxsgf = MAX(maxsgf, imax)
1189 : END IF
1190 : END IF
1191 :
1192 9379097 : IF (PRESENT(maxsgf_set)) THEN
1193 541764 : IF (ASSOCIATED(tmp_basis_set)) THEN
1194 541676 : CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, maxsgf_set=imax)
1195 541676 : maxsgf_set = MAX(maxsgf_set, imax)
1196 : END IF
1197 : END IF
1198 :
1199 9379097 : IF (PRESENT(ncgf)) THEN
1200 48104 : IF (ASSOCIATED(tmp_basis_set)) THEN
1201 13921 : CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, ncgf=n)
1202 13921 : ncgf = ncgf + n*qs_kind_set(ikind)%natom
1203 34183 : ELSE IF (ASSOCIATED(qs_kind%dftb_parameter)) THEN
1204 1602 : CALL get_dftb_atom_param(dftb_parameter=dftb_parameter, lmax=imax)
1205 1602 : n = ((imax + 1)*(imax + 2)*(imax + 3))/6
1206 1602 : ncgf = ncgf + n*qs_kind_set(ikind)%natom
1207 : END IF
1208 : END IF
1209 :
1210 9379097 : IF (PRESENT(npgf)) THEN
1211 42990 : IF (ASSOCIATED(tmp_basis_set)) THEN
1212 8834 : CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, npgf_sum=n)
1213 8834 : npgf = npgf + n*qs_kind_set(ikind)%natom
1214 : END IF
1215 : END IF
1216 :
1217 9379097 : IF (PRESENT(nset)) THEN
1218 42990 : IF (ASSOCIATED(tmp_basis_set)) THEN
1219 8834 : CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, nset=n)
1220 8834 : nset = nset + n*qs_kind_set(ikind)%natom
1221 : END IF
1222 : END IF
1223 :
1224 9379097 : IF (PRESENT(nsgf)) THEN
1225 133905 : IF (ASSOCIATED(tmp_basis_set)) THEN
1226 79406 : CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, nsgf=n)
1227 79406 : nsgf = nsgf + n*qs_kind_set(ikind)%natom
1228 54499 : ELSE IF (ASSOCIATED(qs_kind%dftb_parameter)) THEN
1229 21916 : CALL get_dftb_atom_param(dftb_parameter=dftb_parameter, natorb=n)
1230 21916 : nsgf = nsgf + n*qs_kind_set(ikind)%natom
1231 : END IF
1232 : END IF
1233 :
1234 9379097 : IF (PRESENT(nshell)) THEN
1235 43016 : IF (ASSOCIATED(tmp_basis_set)) THEN
1236 8860 : CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, nshell_sum=n)
1237 8860 : nshell = nshell + n*qs_kind_set(ikind)%natom
1238 34156 : ELSE IF (ASSOCIATED(qs_kind%dftb_parameter)) THEN
1239 1575 : CALL get_dftb_atom_param(dftb_parameter=dftb_parameter, lmax=n)
1240 1575 : nshell = nshell + (n + 1)*qs_kind_set(ikind)%natom
1241 : END IF
1242 : END IF
1243 :
1244 9379097 : IF (PRESENT(nelectron)) THEN
1245 266108 : IF (ASSOCIATED(qs_kind%all_potential)) THEN
1246 : CALL get_potential(potential=qs_kind%all_potential, &
1247 27080 : zeff=zeff, zeff_correction=zeff_correction)
1248 239028 : ELSE IF (ASSOCIATED(qs_kind%gth_potential)) THEN
1249 : CALL get_potential(potential=qs_kind%gth_potential, &
1250 236442 : zeff=zeff, zeff_correction=zeff_correction)
1251 2586 : ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
1252 : CALL get_potential(potential=qs_kind%sgp_potential, &
1253 1610 : zeff=zeff, zeff_correction=zeff_correction)
1254 976 : ELSE IF (ASSOCIATED(qs_kind%cneo_potential)) THEN
1255 : CALL get_cneo_potential(potential=qs_kind%cneo_potential, &
1256 56 : zeff=zeff)
1257 56 : zeff_correction = 0.0_dp
1258 : ELSE
1259 920 : zeff = 0.0_dp
1260 920 : zeff_correction = 0.0_dp
1261 : END IF
1262 266108 : nelectron = nelectron + qs_kind_set(ikind)%natom*NINT(zeff - zeff_correction)
1263 : END IF
1264 :
1265 9379097 : IF (PRESENT(basis_rcut)) THEN
1266 286 : IF (ASSOCIATED(tmp_basis_set)) THEN
1267 0 : CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, kind_radius=brcut)
1268 0 : basis_rcut = MAX(basis_rcut, brcut)
1269 286 : ELSE IF (ASSOCIATED(qs_kind%dftb_parameter)) THEN
1270 286 : CALL get_dftb_atom_param(dftb_parameter=dftb_parameter, cutoff=brcut)
1271 286 : basis_rcut = MAX(basis_rcut, brcut)
1272 : END IF
1273 : END IF
1274 :
1275 9379097 : IF (PRESENT(total_zeff_corr)) THEN
1276 17225 : IF (ASSOCIATED(qs_kind%all_potential)) THEN
1277 : CALL get_potential(potential=qs_kind%all_potential, &
1278 6948 : zeff=zeff, zeff_correction=zeff_correction)
1279 10277 : ELSE IF (ASSOCIATED(qs_kind%gth_potential)) THEN
1280 : CALL get_potential(potential=qs_kind%gth_potential, &
1281 10009 : zeff=zeff, zeff_correction=zeff_correction)
1282 268 : ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
1283 : CALL get_potential(potential=qs_kind%sgp_potential, &
1284 106 : zeff=zeff, zeff_correction=zeff_correction)
1285 : ELSE
1286 162 : zeff = 0.0_dp
1287 162 : zeff_correction = 0.0_dp
1288 : END IF
1289 17225 : total_zeff_corr = total_zeff_corr + qs_kind_set(ikind)%natom*zeff_correction
1290 : END IF
1291 :
1292 9379097 : IF (PRESENT(all_potential_present)) THEN
1293 83591 : IF (ASSOCIATED(all_potential)) THEN
1294 50338 : all_potential_present = .TRUE.
1295 : END IF
1296 : END IF
1297 :
1298 9379097 : IF (PRESENT(tnadd_potential_present)) THEN
1299 0 : IF (ASSOCIATED(tnadd_potential)) THEN
1300 0 : tnadd_potential_present = .TRUE.
1301 : END IF
1302 : END IF
1303 :
1304 9379097 : IF (PRESENT(gth_potential_present)) THEN
1305 66310 : IF (ASSOCIATED(gth_potential)) THEN
1306 22616 : gth_potential_present = .TRUE.
1307 : END IF
1308 : END IF
1309 :
1310 9379097 : IF (PRESENT(sgp_potential_present)) THEN
1311 66316 : IF (ASSOCIATED(sgp_potential)) THEN
1312 164 : sgp_potential_present = .TRUE.
1313 : END IF
1314 : END IF
1315 :
1316 9379097 : IF (PRESENT(cneo_potential_present)) THEN
1317 90497 : IF (ASSOCIATED(cneo_potential)) THEN
1318 24 : cneo_potential_present = .TRUE.
1319 : END IF
1320 : END IF
1321 :
1322 9379097 : IF (PRESENT(nkind_q)) THEN
1323 8598 : IF (ASSOCIATED(cneo_potential)) THEN
1324 4 : nkind_q = nkind_q + 1
1325 : END IF
1326 : END IF
1327 :
1328 9379097 : IF (PRESENT(natom_q)) THEN
1329 8598 : IF (ASSOCIATED(cneo_potential)) THEN
1330 4 : natom_q = natom_q + qs_kind_set(ikind)%natom
1331 : END IF
1332 : END IF
1333 :
1334 9379097 : IF (PRESENT(paw_atom_present)) THEN
1335 66310 : IF (paw_atom) THEN
1336 4902 : paw_atom_present = .TRUE.
1337 : END IF
1338 : END IF
1339 :
1340 9379097 : IF (PRESENT(dft_plus_u_atom_present)) THEN
1341 17225 : IF (dft_plus_u_atom) THEN
1342 90 : dft_plus_u_atom_present = .TRUE.
1343 : END IF
1344 : END IF
1345 :
1346 9379097 : IF (PRESENT(do_mtlr_present)) THEN
1347 17225 : IF (do_mtlr) THEN
1348 16 : do_mtlr_present = .TRUE.
1349 : END IF
1350 : END IF
1351 :
1352 9379097 : IF (PRESENT(max_ngrid_rad)) THEN
1353 0 : max_ngrid_rad = MAX(max_ngrid_rad, ngrid_rad)
1354 : END IF
1355 :
1356 9379097 : IF (PRESENT(max_sph_harm)) THEN
1357 0 : max_sph_harm = MAX(max_sph_harm, max_s_harm)
1358 : END IF
1359 :
1360 9379097 : IF (PRESENT(maxg_iso_not0)) THEN
1361 55358 : maxg_iso_not0 = MAX(maxg_iso_not0, max_iso_not0)
1362 : END IF
1363 :
1364 9379097 : IF (PRESENT(lmax_rho0)) THEN
1365 0 : lmax_rho0 = MAX(lmax_rho0, lmax_rho0_kind)
1366 : END IF
1367 :
1368 23430180 : IF (PRESENT(npgf_seg)) THEN
1369 26 : IF (ASSOCIATED(tmp_basis_set)) THEN
1370 26 : CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, npgf_seg_sum=n)
1371 26 : npgf_seg = npgf_seg + n*qs_kind_set(ikind)%natom
1372 : END IF
1373 : END IF
1374 :
1375 : END DO
1376 : ELSE
1377 0 : CPABORT("The pointer qs_kind_set is not associated")
1378 : END IF
1379 :
1380 4671986 : END SUBROUTINE get_qs_kind_set
1381 :
1382 : ! **************************************************************************************************
1383 : !> \brief Initialise an atomic kind data set.
1384 : !> \param qs_kind ...
1385 : !> \author Creation (11.01.2002,MK)
1386 : !> 20.09.2002 adapted for pol/kg use, gtb
1387 : ! **************************************************************************************************
1388 17225 : SUBROUTINE init_qs_kind(qs_kind)
1389 : TYPE(qs_kind_type), POINTER :: qs_kind
1390 :
1391 : CHARACTER(len=*), PARAMETER :: routineN = 'init_qs_kind'
1392 :
1393 : CHARACTER(LEN=default_string_length) :: basis_type
1394 : INTEGER :: handle, i
1395 : TYPE(gto_basis_set_type), POINTER :: tmp_basis_set
1396 :
1397 17225 : CALL timeset(routineN, handle)
1398 :
1399 17225 : CPASSERT(ASSOCIATED(qs_kind))
1400 :
1401 17225 : IF (ASSOCIATED(qs_kind%gth_potential)) THEN
1402 10009 : CALL init_potential(qs_kind%gth_potential)
1403 7216 : ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
1404 106 : CALL init_potential(qs_kind%sgp_potential)
1405 : END IF
1406 :
1407 361725 : DO i = 1, SIZE(qs_kind%basis_sets, 1)
1408 344500 : NULLIFY (tmp_basis_set)
1409 : CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis_set, &
1410 344500 : inumbas=i, basis_type=basis_type)
1411 344500 : IF (basis_type == "") CYCLE
1412 36691 : IF (basis_type == "AUX") THEN
1413 0 : IF (tmp_basis_set%norm_type < 0) tmp_basis_set%norm_type = 1
1414 0 : CALL init_aux_basis_set(tmp_basis_set)
1415 : ELSE
1416 19466 : IF (tmp_basis_set%norm_type < 0) tmp_basis_set%norm_type = 2
1417 19466 : CALL init_orb_basis_set(tmp_basis_set)
1418 : END IF
1419 : END DO
1420 :
1421 17225 : CALL timestop(handle)
1422 :
1423 17225 : END SUBROUTINE init_qs_kind
1424 :
1425 : ! **************************************************************************************************
1426 : !> \brief Initialise an atomic kind set data set.
1427 : !> \param qs_kind_set ...
1428 : !> \author - Creation (17.01.2002,MK)
1429 : !> - 20.09.2002 para_env passed (gt)
1430 : ! **************************************************************************************************
1431 9324 : SUBROUTINE init_qs_kind_set(qs_kind_set)
1432 :
1433 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1434 :
1435 : CHARACTER(len=*), PARAMETER :: routineN = 'init_qs_kind_set'
1436 :
1437 : INTEGER :: handle, ikind
1438 : TYPE(qs_kind_type), POINTER :: qs_kind
1439 :
1440 9324 : CALL timeset(routineN, handle)
1441 :
1442 9324 : IF (.NOT. ASSOCIATED(qs_kind_set)) THEN
1443 0 : CPABORT("init_qs_kind_set: The pointer qs_kind_set is not associated")
1444 : END IF
1445 :
1446 26549 : DO ikind = 1, SIZE(qs_kind_set)
1447 17225 : qs_kind => qs_kind_set(ikind)
1448 26549 : CALL init_qs_kind(qs_kind)
1449 : END DO
1450 :
1451 9324 : CALL timestop(handle)
1452 :
1453 9324 : END SUBROUTINE init_qs_kind_set
1454 :
1455 : ! **************************************************************************************************
1456 : !> \brief Resolve the active model's pseudopotential GAPW representation.
1457 : !> \param force_env_section ...
1458 : !> \return CP2K default when no GauXC model replaces the XC functional
1459 : ! **************************************************************************************************
1460 1702 : FUNCTION skala_gapw_representation_from_input(force_env_section) RESULT(representation)
1461 : TYPE(section_vals_type), POINTER :: force_env_section
1462 : INTEGER :: representation
1463 :
1464 : CHARACTER(LEN=default_path_length) :: model_key, model_name
1465 : CHARACTER(LEN=default_string_length) :: functional_key, functional_name
1466 : INTEGER :: ifun
1467 : TYPE(section_vals_type), POINTER :: functional, functionals
1468 :
1469 1702 : representation = skala_gapw_cp2k_default
1470 1702 : functionals => section_vals_get_subs_vals(force_env_section, "DFT%XC%XC_FUNCTIONAL")
1471 1702 : IF (.NOT. ASSOCIATED(functionals)) RETURN
1472 :
1473 1702 : ifun = 0
1474 : DO
1475 3424 : ifun = ifun + 1
1476 3424 : functional => section_vals_get_subs_vals2(functionals, i_section=ifun)
1477 3424 : IF (.NOT. ASSOCIATED(functional)) EXIT
1478 1796 : IF (functional%section%name /= "GAUXC") CYCLE
1479 :
1480 74 : CALL section_vals_val_get(functional, "MODEL", c_val=model_name)
1481 74 : CALL section_vals_val_get(functional, "FUNCTIONAL", c_val=functional_name)
1482 74 : model_key = ADJUSTL(model_name)
1483 74 : functional_key = ADJUSTL(functional_name)
1484 74 : CALL uppercase(model_key)
1485 74 : CALL uppercase(functional_key)
1486 74 : IF (TRIM(model_key) /= "" .AND. TRIM(model_key) /= "NONE" .AND. &
1487 : TRIM(model_key) /= TRIM(functional_key)) THEN
1488 : CALL section_vals_val_get(functional, "PSEUDOPOTENTIAL_GAPW_REPRESENTATION", &
1489 66 : i_val=representation)
1490 : END IF
1491 3424 : RETURN
1492 : END DO
1493 :
1494 : END FUNCTION skala_gapw_representation_from_input
1495 :
1496 : ! **************************************************************************************************
1497 : !> \brief ...
1498 : !> \param qs_kind_set ...
1499 : !> \param qs_control ...
1500 : !> \param force_env_section ...
1501 : !> \param modify_qs_control whether the qs_control should be modified
1502 : ! **************************************************************************************************
1503 1702 : SUBROUTINE init_gapw_basis_set(qs_kind_set, qs_control, force_env_section, modify_qs_control)
1504 :
1505 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1506 : TYPE(qs_control_type), POINTER :: qs_control
1507 : TYPE(section_vals_type), POINTER :: force_env_section
1508 : LOGICAL, OPTIONAL :: modify_qs_control
1509 :
1510 : CHARACTER(LEN=default_string_length) :: bsname
1511 : INTEGER :: bas1c, ikind, ilevel, nkind, &
1512 : representation
1513 : LOGICAL :: effective_force_paw, effective_gpw, gpw, &
1514 : has_pseudopotential, my_mod_control, &
1515 : paw_atom
1516 : REAL(dp) :: max_rad_local_type, rc
1517 : TYPE(gth_potential_type), POINTER :: gth_potential
1518 : TYPE(gto_basis_set_type), POINTER :: basis_1c, orb_basis, soft_basis
1519 : TYPE(paw_proj_set_type), POINTER :: paw_proj
1520 : TYPE(qs_kind_type), POINTER :: qs_kind
1521 : TYPE(sgp_potential_type), POINTER :: sgp_potential
1522 :
1523 1702 : my_mod_control = .TRUE.
1524 1702 : IF (PRESENT(modify_qs_control)) THEN
1525 152 : my_mod_control = modify_qs_control
1526 : END IF
1527 1702 : representation = skala_gapw_representation_from_input(force_env_section)
1528 :
1529 1702 : IF (ASSOCIATED(qs_kind_set)) THEN
1530 :
1531 1702 : IF (my_mod_control) qs_control%gapw_control%non_paw_atoms = .FALSE.
1532 1702 : nkind = SIZE(qs_kind_set)
1533 :
1534 4738 : DO ikind = 1, nkind
1535 :
1536 3036 : qs_kind => qs_kind_set(ikind)
1537 :
1538 3036 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis)
1539 3036 : NULLIFY (gth_potential, sgp_potential)
1540 : CALL get_qs_kind(qs_kind=qs_kind, hard_radius=rc, &
1541 : max_rad_local=max_rad_local_type, gpw_type_forced=gpw, &
1542 3036 : gth_potential=gth_potential, sgp_potential=sgp_potential)
1543 :
1544 3036 : has_pseudopotential = ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)
1545 3036 : effective_gpw = gpw
1546 : effective_force_paw = qs_control%gapw_control%force_paw .OR. &
1547 3036 : ASSOCIATED(qs_kind%cneo_potential)
1548 3036 : IF (has_pseudopotential) THEN
1549 20 : SELECT CASE (representation)
1550 : CASE (skala_gapw_direct_valence)
1551 20 : effective_gpw = .TRUE.
1552 20 : effective_force_paw = .FALSE.
1553 : CASE (skala_gapw_paw_one_center, skala_gapw_paw_one_center_split)
1554 44 : effective_gpw = .FALSE.
1555 44 : effective_force_paw = .TRUE.
1556 : CASE (skala_gapw_cp2k_default)
1557 0 : CONTINUE
1558 : CASE DEFAULT
1559 1482 : CPABORT("Unknown GAUXC PSEUDOPOTENTIAL_GAPW_REPRESENTATION value")
1560 : END SELECT
1561 1482 : IF (ASSOCIATED(qs_kind%cneo_potential)) THEN
1562 0 : effective_gpw = .FALSE.
1563 0 : effective_force_paw = .TRUE.
1564 : END IF
1565 1482 : qs_kind%gpw_type_forced = effective_gpw
1566 : END IF
1567 :
1568 3036 : NULLIFY (soft_basis)
1569 3036 : CALL allocate_gto_basis_set(soft_basis)
1570 : ! Quantum nuclear wave functions are very localized. Even if soft
1571 : ! electronic basis is used, the atomic kind needs PAW treatment
1572 : ! because the local Hartree potential is needed.
1573 : CALL create_soft_basis(orb_basis, soft_basis, &
1574 : qs_control%gapw_control%eps_fit, rc, paw_atom, &
1575 3036 : effective_force_paw, effective_gpw)
1576 3036 : CALL add_basis_set_to_container(qs_kind%basis_sets, soft_basis, "ORB_SOFT")
1577 3036 : CALL set_qs_kind(qs_kind=qs_kind, paw_atom=paw_atom)
1578 :
1579 3036 : bas1c = qs_control%gapw_control%basis_1c
1580 3036 : NULLIFY (basis_1c)
1581 2824 : SELECT CASE (bas1c)
1582 : CASE (gapw_1c_orb)
1583 2824 : ilevel = 0
1584 : CASE (gapw_1c_small)
1585 196 : ilevel = 1
1586 : CASE (gapw_1c_medium)
1587 4 : ilevel = 2
1588 : CASE (gapw_1c_large)
1589 8 : ilevel = 3
1590 : CASE (gapw_1c_very_large)
1591 4 : ilevel = 4
1592 : CASE DEFAULT
1593 3036 : CPABORT("basis_1c type")
1594 : END SELECT
1595 3036 : CALL remove_basis_from_container(qs_kind%basis_sets, basis_type="GAPW_1C")
1596 3036 : CALL create_1c_basis(orb_basis, soft_basis, basis_1c, ilevel)
1597 3036 : CALL get_gto_basis_set(gto_basis_set=orb_basis, name=bsname)
1598 3036 : basis_1c%name = TRIM(bsname)//"_1c"
1599 3036 : CALL add_basis_set_to_container(qs_kind%basis_sets, basis_1c, "GAPW_1C")
1600 3036 : IF (paw_atom) THEN
1601 2548 : CALL allocate_paw_proj_set(qs_kind%paw_proj_set)
1602 2548 : CALL get_qs_kind(qs_kind=qs_kind, paw_proj_set=paw_proj)
1603 : CALL projectors(paw_proj, basis_1c, orb_basis, rc, qs_control, &
1604 2548 : max_rad_local_type, force_env_section)
1605 : ELSE
1606 488 : IF (my_mod_control) qs_control%gapw_control%non_paw_atoms = .TRUE.
1607 : END IF
1608 :
1609 : ! grid_atom and harmonics are allocated even if NOT PAW_ATOM
1610 3036 : NULLIFY (qs_kind%grid_atom, qs_kind%harmonics)
1611 3036 : CALL allocate_grid_atom(qs_kind%grid_atom)
1612 10810 : CALL allocate_harmonics_atom(qs_kind%harmonics)
1613 :
1614 : END DO
1615 :
1616 1702 : IF (my_mod_control) THEN
1617 1550 : IF (qs_control%gapw_control%non_paw_atoms) THEN
1618 270 : qs_control%gapw_control%nopaw_as_gpw = .TRUE.
1619 : ELSE
1620 1280 : qs_control%gapw_control%nopaw_as_gpw = .FALSE.
1621 : END IF
1622 : END IF
1623 : ELSE
1624 0 : CPABORT("The pointer qs_kind_set is not associated")
1625 : END IF
1626 :
1627 1702 : END SUBROUTINE init_gapw_basis_set
1628 : ! **************************************************************************************************
1629 : !> \brief ...
1630 : !> \param qs_kind_set ...
1631 : ! **************************************************************************************************
1632 1702 : SUBROUTINE init_gapw_nlcc(qs_kind_set)
1633 :
1634 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1635 :
1636 : INTEGER :: i, ic, ikind, n_nlcc, nc, nexp_nlcc, &
1637 : nkind, nr
1638 1702 : INTEGER, DIMENSION(:), POINTER :: nct_nlcc
1639 : LOGICAL :: nlcc, nlcc_type, paw_atom
1640 : REAL(dp) :: alpha, coa, cval
1641 1702 : REAL(KIND=dp), DIMENSION(:), POINTER :: a_nlcc, alpha_nlcc, c_nlcc, fe, rc, rr
1642 1702 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cval_nlcc, den
1643 : TYPE(gth_potential_type), POINTER :: gth_potential
1644 : TYPE(qs_kind_type), POINTER :: qs_kind
1645 : TYPE(sgp_potential_type), POINTER :: sgp_potential
1646 :
1647 1702 : IF (ASSOCIATED(qs_kind_set)) THEN
1648 1702 : nlcc = has_nlcc(qs_kind_set)
1649 1702 : IF (nlcc) THEN
1650 28 : nkind = SIZE(qs_kind_set)
1651 82 : DO ikind = 1, nkind
1652 54 : qs_kind => qs_kind_set(ikind)
1653 54 : CALL get_qs_kind(qs_kind, paw_atom=paw_atom)
1654 82 : IF (paw_atom) THEN
1655 50 : CALL get_qs_kind(qs_kind, gth_potential=gth_potential)
1656 50 : CALL get_qs_kind(qs_kind, sgp_potential=sgp_potential)
1657 50 : IF (ASSOCIATED(gth_potential)) THEN
1658 : CALL get_potential(potential=gth_potential, nlcc_present=nlcc_type, &
1659 42 : nexp_nlcc=nexp_nlcc, alpha_nlcc=alpha_nlcc, nct_nlcc=nct_nlcc, cval_nlcc=cval_nlcc)
1660 42 : IF (nlcc_type) THEN
1661 24 : nr = qs_kind%grid_atom%nr
1662 24 : rr => qs_kind%grid_atom%rad
1663 144 : ALLOCATE (qs_kind%nlcc_pot(nr, 2), rc(nr), fe(nr))
1664 72 : den => qs_kind%nlcc_pot
1665 2472 : den = 0.0_dp
1666 48 : DO i = 1, nexp_nlcc
1667 24 : alpha = alpha_nlcc(i)
1668 2424 : rc(:) = rr(:)/alpha
1669 2424 : fe(:) = EXP(-0.5_dp*rc(:)*rc(:))
1670 24 : nc = nct_nlcc(i)
1671 74 : DO ic = 1, nc
1672 26 : cval = cval_nlcc(ic, i)
1673 26 : coa = cval/alpha
1674 2626 : den(:, 1) = den(:, 1) + fe(:)*rc**(2*ic - 2)*cval
1675 2626 : den(:, 2) = den(:, 2) - fe(:)*rc**(2*ic - 1)*coa
1676 50 : IF (ic > 1) THEN
1677 202 : den(:, 2) = den(:, 2) + REAL(2*ic - 2, dp)*fe(:)*rc**(2*ic - 3)*coa
1678 : END IF
1679 : END DO
1680 : END DO
1681 24 : DEALLOCATE (rc, fe)
1682 : END IF
1683 8 : ELSE IF (ASSOCIATED(sgp_potential)) THEN
1684 : CALL get_potential(potential=sgp_potential, has_nlcc=nlcc_type, &
1685 8 : n_nlcc=n_nlcc, a_nlcc=a_nlcc, c_nlcc=c_nlcc)
1686 8 : IF (nlcc_type) THEN
1687 4 : nr = qs_kind%grid_atom%nr
1688 4 : rr => qs_kind%grid_atom%rad
1689 24 : ALLOCATE (qs_kind%nlcc_pot(nr, 2), rc(nr), fe(nr))
1690 12 : den => qs_kind%nlcc_pot
1691 652 : den = 0.0_dp
1692 196 : DO i = 1, n_nlcc
1693 192 : alpha = a_nlcc(i)
1694 30912 : fe(:) = EXP(-alpha*rr(:)*rr(:))
1695 192 : cval = c_nlcc(i)
1696 30912 : den(:, 1) = den(:, 1) + cval*fe(:)
1697 30916 : den(:, 2) = den(:, 2) - 2.0_dp*alpha*cval*rr(:)*fe(:)
1698 : END DO
1699 4 : DEALLOCATE (rc, fe)
1700 : END IF
1701 : ELSE
1702 : ! skip
1703 : END IF
1704 : END IF
1705 : END DO
1706 : END IF
1707 : ELSE
1708 0 : CPABORT("The pointer qs_kind_set is not associated")
1709 : END IF
1710 :
1711 1702 : END SUBROUTINE init_gapw_nlcc
1712 :
1713 : ! **************************************************************************************************
1714 : !> \brief ...
1715 : !> \param qs_kind_set ...
1716 : !> \param qs_control ...
1717 : ! **************************************************************************************************
1718 8 : SUBROUTINE init_cneo_basis_set(qs_kind_set, qs_control)
1719 :
1720 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1721 : TYPE(qs_control_type), POINTER :: qs_control
1722 :
1723 : INTEGER :: ikind, nkind
1724 : LOGICAL :: paw_atom
1725 : REAL(dp) :: rc
1726 : TYPE(gto_basis_set_type), POINTER :: orb_basis, soft_basis
1727 : TYPE(qs_kind_type), POINTER :: qs_kind
1728 :
1729 8 : IF (ASSOCIATED(qs_kind_set)) THEN
1730 :
1731 8 : nkind = SIZE(qs_kind_set)
1732 :
1733 22 : DO ikind = 1, nkind
1734 :
1735 14 : qs_kind => qs_kind_set(ikind)
1736 :
1737 22 : IF (ASSOCIATED(qs_kind%cneo_potential)) THEN
1738 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis, basis_type="NUC", &
1739 8 : hard_radius=rc)
1740 :
1741 8 : NULLIFY (soft_basis)
1742 8 : CALL allocate_gto_basis_set(soft_basis)
1743 : CALL create_soft_basis(orb_basis, soft_basis, &
1744 : qs_control%gapw_control%eps_fit/ &
1745 : SQRT(qs_kind%cneo_potential%zeff), &
1746 8 : rc, paw_atom, .TRUE., .FALSE.)
1747 8 : CALL add_basis_set_to_container(qs_kind%basis_sets, soft_basis, "NUC_SOFT")
1748 : END IF
1749 :
1750 : END DO
1751 :
1752 : ELSE
1753 0 : CPABORT("The pointer qs_kind_set is not associated")
1754 : END IF
1755 :
1756 8 : END SUBROUTINE init_cneo_basis_set
1757 :
1758 : ! **************************************************************************************************
1759 : !> \brief Select the DFT potential-file list or a per-kind override.
1760 : !> \param potential_file_names ...
1761 : !> \param dft_potential_file_names ...
1762 : !> \param kind_potential_file_name ...
1763 : ! **************************************************************************************************
1764 34670 : SUBROUTINE set_potential_file_names(potential_file_names, dft_potential_file_names, &
1765 : kind_potential_file_name)
1766 :
1767 : CHARACTER(LEN=default_path_length), DIMENSION(:), &
1768 : INTENT(INOUT), POINTER :: potential_file_names
1769 : CHARACTER(LEN=default_path_length), DIMENSION(:), &
1770 : INTENT(IN) :: dft_potential_file_names
1771 : CHARACTER(LEN=*), INTENT(IN) :: kind_potential_file_name
1772 :
1773 17335 : IF (ASSOCIATED(potential_file_names)) DEALLOCATE (potential_file_names)
1774 17335 : IF (kind_potential_file_name == "-") THEN
1775 51975 : ALLOCATE (potential_file_names(SIZE(dft_potential_file_names)))
1776 34670 : potential_file_names(:) = dft_potential_file_names(:)
1777 : ELSE
1778 10 : ALLOCATE (potential_file_names(1))
1779 10 : potential_file_names(1) = kind_potential_file_name
1780 : END IF
1781 :
1782 17335 : END SUBROUTINE set_potential_file_names
1783 :
1784 : ! **************************************************************************************************
1785 : !> \brief Read an atomic kind data set from the input file.
1786 : !> \param qs_kind ...
1787 : !> \param kind_section ...
1788 : !> \param para_env ...
1789 : !> \param force_env_section ...
1790 : !> \param no_fail ...
1791 : !> \param method_id ...
1792 : !> \param silent ...
1793 : !> \par History
1794 : !> - Creation (09.02.2002,MK)
1795 : !> - 20.09.2002,gt: adapted for POL/KG use (elp_potential)
1796 : !> - 05.03.2010: split elp_potential into fist_potential and kg_potential
1797 : ! **************************************************************************************************
1798 17311 : SUBROUTINE read_qs_kind(qs_kind, kind_section, para_env, force_env_section, &
1799 : no_fail, method_id, silent)
1800 :
1801 : TYPE(qs_kind_type), INTENT(INOUT) :: qs_kind
1802 : TYPE(section_vals_type), POINTER :: kind_section
1803 : TYPE(mp_para_env_type), POINTER :: para_env
1804 : TYPE(section_vals_type), POINTER :: force_env_section
1805 : LOGICAL, INTENT(IN) :: no_fail
1806 : INTEGER, INTENT(IN) :: method_id
1807 : LOGICAL, INTENT(IN) :: silent
1808 :
1809 : CHARACTER(LEN=*), PARAMETER :: routineN = 'read_qs_kind'
1810 : INTEGER, PARAMETER :: maxbas = 20
1811 :
1812 : CHARACTER(LEN=2) :: element_symbol
1813 : CHARACTER(len=default_path_length) :: kg_potential_fn_kind, potential_fn_kind
1814 : CHARACTER(LEN=default_path_length), DIMENSION(:), &
1815 17311 : POINTER :: dft_potential_file_names, &
1816 17311 : potential_file_names
1817 : CHARACTER(LEN=default_string_length) :: akind_name, basis_type, keyword, &
1818 : kgpot_name, kgpot_type, &
1819 : potential_name, potential_type, tmp
1820 : CHARACTER(LEN=default_string_length), DIMENSION(4) :: description
1821 : CHARACTER(LEN=default_string_length), &
1822 17311 : DIMENSION(:), POINTER :: tmpstringlist
1823 : CHARACTER(LEN=default_string_length), &
1824 : DIMENSION(maxbas) :: basis_set_form, basis_set_name, &
1825 : basis_set_type
1826 : INTEGER :: handle, i, i_rep, iounit, ipaodesc, ipaopot, ipos, j, jj, k_rep, l, m, &
1827 : n_potential_files, n_rep, nb_rep, nexp, ngauss, nlcc, nloc, nnl, norbitals, npaodesc, &
1828 : npaopot, nppnl, nspin, nu, plus_u_method_id, z
1829 34622 : INTEGER, DIMENSION(:), POINTER :: add_el, elec_conf, orbitals
1830 : LOGICAL :: check, ecp_semi_local, explicit, explicit_basis, explicit_J, explicit_kgpot, &
1831 : explicit_potential, explicit_U, explicit_u_m_j, l_explicit, n_explicit, nobasis, &
1832 : nobasis_nuc, section_enabled, subsection_enabled, update_input
1833 : REAL(KIND=dp) :: alpha, ccore, mass, r, rc, &
1834 : zeff_correction
1835 : REAL(KIND=dp), DIMENSION(6) :: error
1836 34622 : REAL(KIND=dp), DIMENSION(:), POINTER :: a_nl, aloc, anlcc, cloc, cnlcc, nelec, &
1837 17311 : proj_shell_charge, &
1838 17311 : tmp_perturbation_strength
1839 17311 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: h_nl
1840 17311 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: c_nl
1841 : TYPE(atom_ecppot_type) :: ecppot
1842 : TYPE(atom_sgp_potential_type) :: sgppot
1843 1817655 : TYPE(atom_upfpot_type) :: upfpot
1844 : TYPE(cp_logger_type), POINTER :: logger
1845 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set, sup_basis_set, &
1846 : tmp_basis_set
1847 : TYPE(section_vals_type), POINTER :: basis_section, bs_section, dft_plus_u_section, &
1848 : dft_section, enforce_occupation_section, kgpot_section, &
1849 : minimum_tracking_linear_response_section, pao_desc_section, pao_pot_section, &
1850 : potential_section, spin_section, tensorial_projector_section
1851 : TYPE(sto_basis_set_type), POINTER :: sto_basis_set
1852 :
1853 17311 : CALL timeset(routineN, handle)
1854 :
1855 17311 : NULLIFY (logger)
1856 17311 : logger => cp_get_default_logger()
1857 17311 : iounit = cp_logger_get_default_io_unit(logger)
1858 :
1859 17311 : NULLIFY (dft_potential_file_names, elec_conf, potential_file_names)
1860 :
1861 17311 : update_input = .TRUE.
1862 363531 : basis_set_name(:) = ""
1863 363531 : basis_set_type(:) = ""
1864 363531 : basis_set_form(:) = ""
1865 17311 : potential_name = ""
1866 17311 : potential_type = ""
1867 17311 : kgpot_name = ""
1868 17311 : kgpot_type = ""
1869 17311 : kg_potential_fn_kind = "-"
1870 17311 : potential_fn_kind = "-"
1871 17311 : z = -1
1872 17311 : zeff_correction = 0.0_dp
1873 17311 : explicit = .FALSE.
1874 17311 : explicit_basis = .FALSE.
1875 17311 : explicit_J = .FALSE.
1876 17311 : explicit_kgpot = .FALSE.
1877 17311 : explicit_potential = .FALSE.
1878 17311 : explicit_U = .FALSE.
1879 17311 : explicit_u_m_j = .FALSE.
1880 :
1881 17311 : dft_section => section_vals_get_subs_vals(force_env_section, "DFT")
1882 17311 : CALL section_vals_val_get(dft_section, "POTENTIAL_FILE_NAME", n_rep_val=n_potential_files)
1883 17311 : CPASSERT(n_potential_files > 0)
1884 51933 : ALLOCATE (dft_potential_file_names(n_potential_files))
1885 34642 : DO i = 1, n_potential_files
1886 : CALL section_vals_val_get(dft_section, "POTENTIAL_FILE_NAME", i_rep_val=i, &
1887 34642 : c_val=dft_potential_file_names(i))
1888 : END DO
1889 17311 : CALL section_vals_get(kind_section, n_repetition=n_rep)
1890 17311 : k_rep = -1
1891 17311 : akind_name = qs_kind%name
1892 17311 : CALL uppercase(akind_name)
1893 : ! First we use the atom_name to find out the proper KIND section
1894 24692 : DO i_rep = 1, n_rep
1895 : CALL section_vals_val_get(kind_section, "_SECTION_PARAMETERS_", &
1896 19290 : c_val=keyword, i_rep_section=i_rep)
1897 19290 : CALL uppercase(keyword)
1898 24692 : IF (keyword == akind_name) THEN
1899 11909 : k_rep = i_rep
1900 11909 : EXIT
1901 : END IF
1902 : END DO
1903 : ! The search for the KIND section failed.. check for a QM/MM link atom
1904 17311 : IF (k_rep < 1) THEN
1905 5402 : ipos = INDEX(qs_kind%name, "_")
1906 5402 : IF (((ipos == 2) .OR. (ipos == 3)) .AND. (INDEX(qs_kind%name, "_ghost") == 0)) THEN
1907 : ! If the atm_name could not match any KIND section it maybe be a QM/MM link atom.
1908 : ! ghost atoms will be treated differently.
1909 64 : akind_name = qs_kind%name(1:ipos - 1)
1910 64 : CALL uppercase(akind_name)
1911 64 : DO i_rep = 1, n_rep
1912 : CALL section_vals_val_get(kind_section, "_SECTION_PARAMETERS_", &
1913 52 : c_val=keyword, i_rep_section=i_rep)
1914 52 : CALL uppercase(keyword)
1915 64 : IF (keyword == akind_name) THEN
1916 52 : k_rep = i_rep
1917 52 : EXIT
1918 : END IF
1919 : END DO
1920 : END IF
1921 : END IF
1922 : ! The search for the KIND section failed.. check element_symbol
1923 17311 : IF (k_rep < 1) THEN
1924 : ! If it's not a link atom let's check for the element and map
1925 : ! the KIND section to the element.
1926 5350 : element_symbol = qs_kind%element_symbol(1:2)
1927 5350 : CALL uppercase(element_symbol)
1928 5448 : DO i_rep = 1, n_rep
1929 : CALL section_vals_val_get(kind_section, "_SECTION_PARAMETERS_", &
1930 122 : c_val=keyword, i_rep_section=i_rep)
1931 122 : CALL uppercase(keyword)
1932 5448 : IF (keyword == element_symbol) THEN
1933 24 : k_rep = i_rep
1934 24 : EXIT
1935 : END IF
1936 : END DO
1937 : END IF
1938 : ! In case it should not really match any possible KIND section
1939 : ! let's look if a default one is defined..
1940 17311 : IF (k_rep < 1) THEN
1941 5342 : DO i_rep = 1, n_rep
1942 : CALL section_vals_val_get(kind_section, "_SECTION_PARAMETERS_", &
1943 76 : c_val=keyword, i_rep_section=i_rep)
1944 76 : CALL uppercase(keyword)
1945 5342 : IF (keyword == "DEFAULT") THEN
1946 60 : update_input = .FALSE.
1947 60 : k_rep = i_rep
1948 60 : EXIT
1949 : END IF
1950 : END DO
1951 : END IF
1952 17311 : IF (k_rep < 0 .AND. (.NOT. no_fail)) THEN
1953 : CALL cp_abort(__LOCATION__, &
1954 : "No &KIND section was possible to associate to the atomic kind <"// &
1955 : TRIM(akind_name)//">. The KIND section were also scanned for the"// &
1956 : " corresponding element <"//TRIM(qs_kind%element_symbol)//">"// &
1957 0 : " and for the DEFAULT section but no match was found. Check your input file!")
1958 : END IF
1959 : ! Retrieve information on element
1960 17311 : CALL get_ptable_info(qs_kind%element_symbol, ielement=z)
1961 :
1962 : ! Normal parsing of the KIND section
1963 17311 : IF (k_rep > 0) THEN
1964 : ! new style basis set input
1965 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1966 : keyword_name="BASIS_SET", &
1967 : explicit=explicit, &
1968 12045 : n_rep_val=nb_rep)
1969 12045 : IF (.NOT. explicit) nb_rep = 0
1970 12045 : CPASSERT(nb_rep <= maxbas)
1971 26293 : DO i = 1, nb_rep
1972 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1973 14248 : keyword_name="BASIS_SET", i_rep_val=i, c_vals=tmpstringlist)
1974 14248 : IF (SIZE(tmpstringlist) == 1) THEN
1975 : ! default is orbital type and GTO
1976 10555 : basis_set_type(i) = "ORB"
1977 10555 : basis_set_form(i) = "GTO"
1978 10555 : basis_set_name(i) = tmpstringlist(1)
1979 3693 : ELSE IF (SIZE(tmpstringlist) == 2) THEN
1980 : ! default is GTO
1981 3689 : basis_set_type(i) = tmpstringlist(1)
1982 3689 : basis_set_form(i) = "GTO"
1983 3689 : basis_set_name(i) = tmpstringlist(2)
1984 4 : ELSE IF (SIZE(tmpstringlist) == 3) THEN
1985 4 : basis_set_type(i) = tmpstringlist(1)
1986 4 : basis_set_form(i) = tmpstringlist(2)
1987 4 : basis_set_name(i) = tmpstringlist(3)
1988 : ELSE
1989 : CALL cp_abort(__LOCATION__, &
1990 0 : "invalid number of BASIS_SET keyword parameters: BASIS_SET [<TYPE>] [<FORM>] <NAME>")
1991 : END IF
1992 : ! check that we have a valid basis set form
1993 26293 : IF (basis_set_form(i) /= "GTO" .AND. basis_set_form(i) /= "STO") THEN
1994 0 : CPABORT("invalid BASIS_SET FORM parameter")
1995 : END IF
1996 : END DO
1997 :
1998 : ! parse PAO keywords
1999 : CALL section_vals_val_get(kind_section, keyword_name="PAO_BASIS_SIZE", i_rep_section=k_rep, &
2000 12045 : i_val=qs_kind%pao_basis_size)
2001 : CALL section_vals_val_get(kind_section, keyword_name="PAO_MODEL_FILE", i_rep_section=k_rep, &
2002 12045 : explicit=explicit)
2003 12045 : IF (explicit) THEN
2004 : CALL section_vals_val_get(kind_section, keyword_name="PAO_MODEL_FILE", i_rep_section=k_rep, &
2005 8 : c_val=qs_kind%pao_model_file)
2006 : END IF
2007 :
2008 : ! parse PAO_POTENTIAL sections
2009 12045 : pao_pot_section => section_vals_get_subs_vals(kind_section, "PAO_POTENTIAL", i_rep_section=k_rep)
2010 12045 : CALL section_vals_get(pao_pot_section, n_repetition=npaopot)
2011 24172 : ALLOCATE (qs_kind%pao_potentials(npaopot))
2012 12087 : DO ipaopot = 1, npaopot
2013 : CALL section_vals_val_get(pao_pot_section, keyword_name="MAXL", i_rep_section=ipaopot, &
2014 42 : i_val=qs_kind%pao_potentials(ipaopot)%maxl)
2015 : CALL section_vals_val_get(pao_pot_section, keyword_name="MAX_PROJECTOR", i_rep_section=ipaopot, &
2016 42 : i_val=qs_kind%pao_potentials(ipaopot)%max_projector)
2017 : CALL section_vals_val_get(pao_pot_section, keyword_name="BETA", i_rep_section=ipaopot, &
2018 42 : r_val=qs_kind%pao_potentials(ipaopot)%beta)
2019 : CALL section_vals_val_get(pao_pot_section, keyword_name="WEIGHT", i_rep_section=ipaopot, &
2020 12087 : r_val=qs_kind%pao_potentials(ipaopot)%weight)
2021 : END DO
2022 :
2023 : ! parse PAO_DESCRIPTOR sections
2024 12045 : pao_desc_section => section_vals_get_subs_vals(kind_section, "PAO_DESCRIPTOR", i_rep_section=k_rep)
2025 12045 : CALL section_vals_get(pao_desc_section, n_repetition=npaodesc)
2026 24120 : ALLOCATE (qs_kind%pao_descriptors(npaodesc))
2027 12063 : DO ipaodesc = 1, npaodesc
2028 : CALL section_vals_val_get(pao_desc_section, keyword_name="BETA", i_rep_section=ipaodesc, &
2029 18 : r_val=qs_kind%pao_descriptors(ipaodesc)%beta)
2030 : CALL section_vals_val_get(pao_desc_section, keyword_name="SCREENING", i_rep_section=ipaodesc, &
2031 18 : r_val=qs_kind%pao_descriptors(ipaodesc)%screening)
2032 : CALL section_vals_val_get(pao_desc_section, keyword_name="WEIGHT", i_rep_section=ipaodesc, &
2033 12063 : r_val=qs_kind%pao_descriptors(ipaodesc)%weight)
2034 : END DO
2035 :
2036 : ! parse ELEC_CONF
2037 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
2038 12045 : keyword_name="ELEC_CONF", n_rep_val=i)
2039 12045 : IF (i > 0) THEN
2040 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
2041 0 : keyword_name="ELEC_CONF", i_vals=elec_conf)
2042 0 : CALL set_qs_kind(qs_kind, elec_conf=elec_conf)
2043 : END IF
2044 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
2045 12045 : keyword_name="CORE_CORRECTION", r_val=zeff_correction)
2046 : ! parse POTENTIAL
2047 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
2048 12045 : keyword_name="POTENTIAL_FILE_NAME", c_val=potential_fn_kind)
2049 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
2050 12045 : keyword_name="POTENTIAL_TYPE", c_val=potential_type)
2051 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
2052 12045 : explicit=explicit, keyword_name="POTENTIAL", c_vals=tmpstringlist)
2053 12045 : IF (explicit) THEN
2054 11721 : IF (SIZE(tmpstringlist) == 1) THEN
2055 : ! old type of input: start of name defines type
2056 11593 : potential_name = tmpstringlist(1)
2057 11593 : IF (potential_type == "") THEN
2058 11593 : ipos = INDEX(potential_name, "-")
2059 11593 : IF (ipos > 1) THEN
2060 10111 : potential_type = potential_name(:ipos - 1)
2061 : ELSE
2062 1482 : potential_type = potential_name
2063 : END IF
2064 : END IF
2065 128 : ELSE IF (SIZE(tmpstringlist) == 2) THEN
2066 128 : potential_type = tmpstringlist(1)
2067 128 : potential_name = tmpstringlist(2)
2068 : ELSE
2069 0 : CPABORT("POTENTIAL input list is not correct")
2070 : END IF
2071 : END IF
2072 12045 : CALL uppercase(potential_type)
2073 :
2074 : ! Parse KG POTENTIAL
2075 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
2076 12045 : keyword_name="KG_POTENTIAL_FILE_NAME", c_val=kg_potential_fn_kind)
2077 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
2078 12045 : keyword_name="KG_POTENTIAL", c_val=kgpot_name)
2079 :
2080 : ! Semi-local vs. full nonlocal form of ECPs
2081 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
2082 12045 : keyword_name="ECP_SEMI_LOCAL", l_val=ecp_semi_local)
2083 :
2084 : ! Assign atomic covalent radius
2085 12045 : qs_kind%covalent_radius = ptable(z)%covalent_radius*bohr
2086 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
2087 12045 : keyword_name="COVALENT_RADIUS", r_val=r)
2088 12045 : IF (r > 0.0_dp) qs_kind%covalent_radius = r
2089 :
2090 : ! Assign atomic van der Waals radius
2091 12045 : qs_kind%vdw_radius = ptable(z)%vdw_radius*bohr
2092 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
2093 12045 : keyword_name="VDW_RADIUS", r_val=r)
2094 12045 : IF (r > 0.0_dp) qs_kind%vdw_radius = r
2095 :
2096 : ! Assign atom dependent defaults, only H special case
2097 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, n_rep_val=i, &
2098 12045 : keyword_name="HARD_EXP_RADIUS")
2099 12045 : IF (i == 0) THEN
2100 11923 : IF (z == 1) THEN
2101 5180 : qs_kind%hard_radius = 1.2_dp
2102 : ELSE
2103 6743 : qs_kind%hard_radius = 0.8_dp*bohr
2104 : END IF
2105 : ELSE
2106 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
2107 122 : keyword_name="HARD_EXP_RADIUS", r_val=qs_kind%hard_radius)
2108 : END IF
2109 :
2110 : ! assign atom dependent defaults, only H special case
2111 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, n_rep_val=i, &
2112 12045 : keyword_name="RHO0_EXP_RADIUS")
2113 12045 : IF (i == 0) THEN
2114 12045 : qs_kind%hard0_radius = qs_kind%hard_radius
2115 : ELSE
2116 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
2117 0 : keyword_name="RHO0_EXP_RADIUS", r_val=qs_kind%hard0_radius)
2118 : END IF
2119 12045 : IF (qs_kind%hard_radius < qs_kind%hard0_radius) THEN
2120 0 : CPABORT("rc0 should be <= rc")
2121 : END IF
2122 :
2123 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
2124 12045 : keyword_name="MAX_RAD_LOCAL", r_val=qs_kind%max_rad_local)
2125 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
2126 12045 : keyword_name="LEBEDEV_GRID", i_val=qs_kind%ngrid_ang)
2127 12045 : IF (qs_kind%ngrid_ang <= 0) THEN
2128 0 : CPABORT("# point lebedev grid < 0")
2129 : END IF
2130 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
2131 12045 : keyword_name="RADIAL_GRID", i_val=qs_kind%ngrid_rad)
2132 12045 : IF (qs_kind%ngrid_rad <= 0) THEN
2133 0 : CPABORT("# point radial grid < 0")
2134 : END IF
2135 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
2136 12045 : keyword_name="GPW_TYPE", l_val=qs_kind%gpw_type_forced)
2137 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
2138 12045 : keyword_name="GHOST", l_val=qs_kind%ghost)
2139 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
2140 12045 : keyword_name="FLOATING_BASIS_CENTER", l_val=qs_kind%floating)
2141 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
2142 12045 : keyword_name="NO_OPTIMIZE", l_val=qs_kind%no_optimize)
2143 :
2144 : ! Magnetization
2145 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
2146 12045 : keyword_name="MAGNETIZATION", r_val=qs_kind%magnetization)
2147 : ! DFTB3 param
2148 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
2149 12045 : keyword_name="DFTB3_PARAM", r_val=qs_kind%dudq_dftb3)
2150 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
2151 12045 : keyword_name="LMAX_DFTB", i_val=qs_kind%lmax_dftb)
2152 :
2153 : ! MAOS
2154 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
2155 12045 : keyword_name="MAO", i_val=qs_kind%mao)
2156 :
2157 : ! Use monovalent pseudopotential
2158 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
2159 12045 : keyword_name="MONOVALENT", l_val=qs_kind%monovalent)
2160 12045 : IF (qs_kind%monovalent .AND. TRIM(potential_type) /= 'GTH') THEN
2161 0 : CPABORT("Monovalent pseudopotentials currently implemented for GTH only!")
2162 : END IF
2163 :
2164 : ! Read the BS subsection of the current atomic kind, if enabled
2165 12045 : NULLIFY (bs_section)
2166 : bs_section => section_vals_get_subs_vals(kind_section, "BS", &
2167 12045 : i_rep_section=k_rep)
2168 : section_enabled = .FALSE.
2169 : CALL section_vals_val_get(bs_section, "_SECTION_PARAMETERS_", &
2170 12045 : l_val=section_enabled)
2171 12045 : IF (section_enabled) THEN
2172 : ! test for conflict with magnetization
2173 66 : IF (qs_kind%magnetization /= 0.0_dp) THEN
2174 : CALL cp_abort(__LOCATION__, "BS Section is in conflict with non-zero magnetization "// &
2175 0 : "for this atom kind.")
2176 : END IF
2177 66 : qs_kind%bs_occupation = .TRUE.
2178 : !Alpha spin
2179 66 : NULLIFY (spin_section)
2180 66 : spin_section => section_vals_get_subs_vals(bs_section, "ALPHA")
2181 66 : CALL section_vals_get(spin_section, explicit=explicit)
2182 66 : IF (explicit) THEN
2183 66 : NULLIFY (add_el)
2184 : CALL section_vals_val_get(spin_section, &
2185 66 : keyword_name="NEL", i_vals=add_el)
2186 66 : CPASSERT(ASSOCIATED(add_el))
2187 198 : ALLOCATE (qs_kind%addel(SIZE(add_el), 2))
2188 358 : qs_kind%addel = 0
2189 146 : qs_kind%addel(1:SIZE(add_el), 1) = add_el(1:SIZE(add_el))
2190 66 : NULLIFY (add_el)
2191 : CALL section_vals_val_get(spin_section, &
2192 66 : keyword_name="L", i_vals=add_el)
2193 66 : CPASSERT(ASSOCIATED(add_el))
2194 66 : CPASSERT(SIZE(add_el) == SIZE(qs_kind%addel, 1))
2195 198 : ALLOCATE (qs_kind%laddel(SIZE(add_el), 2))
2196 358 : qs_kind%laddel = 0
2197 146 : qs_kind%laddel(1:SIZE(add_el), 1) = add_el(1:SIZE(add_el))
2198 198 : ALLOCATE (qs_kind%naddel(SIZE(add_el), 2))
2199 358 : qs_kind%naddel = 0
2200 66 : NULLIFY (add_el)
2201 : CALL section_vals_val_get(spin_section, &
2202 66 : keyword_name="N", n_rep_val=i)
2203 66 : IF (i > 0) THEN
2204 : CALL section_vals_val_get(spin_section, &
2205 66 : keyword_name="N", i_vals=add_el)
2206 66 : IF (SIZE(add_el) == SIZE(qs_kind%addel, 1)) THEN
2207 146 : qs_kind%naddel(1:SIZE(add_el), 1) = add_el(1:SIZE(add_el))
2208 : END IF
2209 : END IF
2210 : END IF
2211 : ! Beta spin
2212 66 : NULLIFY (spin_section)
2213 66 : spin_section => section_vals_get_subs_vals(bs_section, "BETA")
2214 66 : CALL section_vals_get(spin_section, explicit=explicit)
2215 66 : IF (explicit) THEN
2216 66 : NULLIFY (add_el)
2217 : CALL section_vals_val_get(spin_section, &
2218 66 : keyword_name="NEL", i_vals=add_el)
2219 66 : CPASSERT(SIZE(add_el) == SIZE(qs_kind%addel, 1))
2220 146 : qs_kind%addel(1:SIZE(add_el), 2) = add_el(1:SIZE(add_el))
2221 358 : qs_kind%addel(:, :) = qs_kind%addel(:, :)
2222 66 : NULLIFY (add_el)
2223 : CALL section_vals_val_get(spin_section, &
2224 66 : keyword_name="L", i_vals=add_el)
2225 66 : CPASSERT(SIZE(add_el) == SIZE(qs_kind%addel, 1))
2226 146 : qs_kind%laddel(1:SIZE(add_el), 2) = add_el(1:SIZE(add_el))
2227 :
2228 : CALL section_vals_val_get(spin_section, &
2229 66 : keyword_name="N", n_rep_val=i)
2230 66 : IF (i > 0) THEN
2231 66 : NULLIFY (add_el)
2232 : CALL section_vals_val_get(spin_section, &
2233 66 : keyword_name="N", i_vals=add_el)
2234 66 : IF (SIZE(add_el) == SIZE(qs_kind%addel, 1)) THEN
2235 146 : qs_kind%naddel(1:SIZE(add_el), 2) = add_el(1:SIZE(add_el))
2236 : END IF
2237 : END IF
2238 : END IF
2239 : END IF
2240 :
2241 : ! Read the DFT+U subsection of the current atomic kind, if enabled
2242 :
2243 12045 : NULLIFY (dft_plus_u_section)
2244 : dft_plus_u_section => section_vals_get_subs_vals(kind_section, &
2245 : subsection_name="DFT_PLUS_U", &
2246 12045 : i_rep_section=k_rep)
2247 : section_enabled = .FALSE.
2248 : CALL section_vals_val_get(dft_plus_u_section, &
2249 : keyword_name="_SECTION_PARAMETERS_", &
2250 12045 : l_val=section_enabled)
2251 156585 : IF (section_enabled) THEN
2252 450 : ALLOCATE (qs_kind%dft_plus_u)
2253 90 : NULLIFY (qs_kind%dft_plus_u%nelec)
2254 90 : NULLIFY (qs_kind%dft_plus_u%orbitals)
2255 : CALL section_vals_val_get(dft_plus_u_section, &
2256 : keyword_name="L", &
2257 90 : i_val=l)
2258 90 : qs_kind%dft_plus_u%l = l
2259 : #if defined(__SIRIUS)
2260 : CALL section_vals_val_get(dft_plus_u_section, &
2261 : keyword_name="N", &
2262 90 : i_val=nu)
2263 90 : qs_kind%dft_plus_u%n = nu
2264 :
2265 : CALL section_vals_val_get(dft_plus_u_section, &
2266 : keyword_name="U", &
2267 : r_val=qs_kind%dft_plus_u%U, &
2268 90 : explicit=explicit_U)
2269 :
2270 : CALL section_vals_val_get(dft_plus_u_section, &
2271 : keyword_name="J", &
2272 : r_val=qs_kind%dft_plus_u%J, &
2273 90 : explicit=explicit_J)
2274 :
2275 : CALL section_vals_val_get(dft_plus_u_section, &
2276 : keyword_name="alpha", &
2277 90 : r_val=qs_kind%dft_plus_u%alpha)
2278 :
2279 : CALL section_vals_val_get(dft_plus_u_section, &
2280 : keyword_name="beta", &
2281 90 : r_val=qs_kind%dft_plus_u%beta)
2282 :
2283 : CALL section_vals_val_get(dft_plus_u_section, &
2284 : keyword_name="J0", &
2285 90 : r_val=qs_kind%dft_plus_u%J0)
2286 :
2287 : CALL section_vals_val_get(dft_plus_u_section, &
2288 : keyword_name="occupation", &
2289 90 : r_val=qs_kind%dft_plus_u%occupation)
2290 : #else
2291 : nu = 0
2292 : #endif
2293 :
2294 : CALL section_vals_val_get(dft_plus_u_section, &
2295 : keyword_name="U_MINUS_J", &
2296 : r_val=qs_kind%dft_plus_u%u_minus_j_target, &
2297 90 : explicit=explicit_u_m_j)
2298 :
2299 : CALL section_vals_val_get(dft_section, &
2300 : keyword_name="PLUS_U_METHOD", &
2301 90 : i_val=plus_u_method_id)
2302 :
2303 90 : IF (plus_u_method_id == plus_u_tensorial .OR. plus_u_method_id == plus_u_lowdin) THEN
2304 42 : IF (explicit_U) THEN
2305 : CALL cp_abort(__LOCATION__, "Lowdin and tensorial DFT+U+J do not allow U "// &
2306 0 : "to be specified. Use U_MINUS_J and J.")
2307 : END IF
2308 : ELSE
2309 48 : IF ((explicit_U .OR. explicit_J) .AND. explicit_u_m_j) THEN
2310 0 : CPABORT("DFT+U| specifying U or J and U_MINUS_J parameters are mutually exclusive.")
2311 : END IF
2312 : END IF
2313 :
2314 : ! Tensorial DFT+U+J needs N, L, J, and U_MINUS_J. Here, N, L, and J are read.
2315 90 : l_explicit = .FALSE.
2316 90 : n_explicit = .FALSE.
2317 90 : explicit_J = .FALSE.
2318 :
2319 90 : IF (plus_u_method_id == plus_u_tensorial) THEN
2320 : CALL section_vals_val_get(dft_plus_u_section, &
2321 : keyword_name="L", &
2322 : i_val=qs_kind%dft_plus_u%l, &
2323 8 : explicit=l_explicit)
2324 : CALL section_vals_val_get(dft_plus_u_section, &
2325 : keyword_name="N", &
2326 : i_val=qs_kind%dft_plus_u%n, &
2327 8 : explicit=n_explicit)
2328 : CALL section_vals_val_get(dft_plus_u_section, &
2329 : keyword_name="J", &
2330 : r_val=qs_kind%dft_plus_u%J, &
2331 8 : explicit=explicit_J)
2332 8 : qs_kind%dft_plus_u%hund_j = qs_kind%dft_plus_u%J
2333 8 : IF (.NOT. l_explicit .OR. .NOT. n_explicit) THEN
2334 : CALL cp_abort(__LOCATION__, "DFT PLUS_U_METHOD TENSORIAL requires explicit L "// &
2335 0 : "and N keywords in every active DFT_PLUS_U section.")
2336 : END IF
2337 : END IF
2338 :
2339 90 : IF (plus_u_method_id == plus_u_lowdin) THEN
2340 : CALL section_vals_val_get(dft_plus_u_section, &
2341 : keyword_name="J", &
2342 : r_val=qs_kind%dft_plus_u%J, &
2343 34 : explicit=explicit_J)
2344 34 : qs_kind%dft_plus_u%hund_j = qs_kind%dft_plus_u%J
2345 : END IF
2346 :
2347 : CALL section_vals_val_get(dft_plus_u_section, &
2348 : keyword_name="U_RAMPING", &
2349 90 : r_val=qs_kind%dft_plus_u%u_ramping)
2350 : CALL section_vals_val_get(dft_plus_u_section, &
2351 : keyword_name="INIT_U_RAMPING_EACH_SCF", &
2352 90 : l_val=qs_kind%dft_plus_u%init_u_ramping_each_scf)
2353 90 : IF (qs_kind%dft_plus_u%u_ramping > 0.0_dp) THEN
2354 8 : qs_kind%dft_plus_u%u_minus_j = 0.0_dp
2355 : ELSE
2356 82 : qs_kind%dft_plus_u%u_minus_j = qs_kind%dft_plus_u%u_minus_j_target
2357 : END IF
2358 : CALL section_vals_val_get(dft_plus_u_section, &
2359 : keyword_name="EPS_U_RAMPING", &
2360 90 : r_val=qs_kind%dft_plus_u%eps_u_ramping)
2361 :
2362 90 : NULLIFY (enforce_occupation_section)
2363 : enforce_occupation_section => section_vals_get_subs_vals(dft_plus_u_section, &
2364 90 : subsection_name="ENFORCE_OCCUPATION")
2365 : subsection_enabled = .FALSE.
2366 : CALL section_vals_val_get(enforce_occupation_section, &
2367 : keyword_name="_SECTION_PARAMETERS_", &
2368 90 : l_val=subsection_enabled)
2369 90 : IF (subsection_enabled) THEN
2370 4 : NULLIFY (nelec)
2371 : CALL section_vals_val_get(enforce_occupation_section, &
2372 : keyword_name="NELEC", &
2373 4 : r_vals=nelec)
2374 4 : nspin = SIZE(nelec)
2375 12 : ALLOCATE (qs_kind%dft_plus_u%nelec(nspin))
2376 8 : qs_kind%dft_plus_u%nelec(:) = nelec(:)
2377 4 : NULLIFY (orbitals)
2378 : CALL section_vals_val_get(enforce_occupation_section, &
2379 : keyword_name="ORBITALS", &
2380 4 : i_vals=orbitals)
2381 4 : norbitals = SIZE(orbitals)
2382 4 : IF (norbitals <= 0 .OR. norbitals > 2*l + 1) THEN
2383 : CALL cp_abort(__LOCATION__, "DFT+U| Invalid number of ORBITALS specified: "// &
2384 0 : "1 to 2*L+1 integer numbers are expected")
2385 : END IF
2386 12 : ALLOCATE (qs_kind%dft_plus_u%orbitals(norbitals))
2387 16 : qs_kind%dft_plus_u%orbitals(:) = orbitals(:)
2388 4 : NULLIFY (orbitals)
2389 16 : DO m = 1, norbitals
2390 12 : IF (qs_kind%dft_plus_u%orbitals(m) > l) THEN
2391 0 : CPABORT("DFT+U| Invalid orbital magnetic quantum number specified: m > l")
2392 : END IF
2393 12 : IF (qs_kind%dft_plus_u%orbitals(m) < -l) THEN
2394 0 : CPABORT("DFT+U| Invalid orbital magnetic quantum number specified: m < -l")
2395 : END IF
2396 52 : DO j = 1, norbitals
2397 48 : IF (j /= m) THEN
2398 24 : IF (qs_kind%dft_plus_u%orbitals(j) == qs_kind%dft_plus_u%orbitals(m)) THEN
2399 0 : CPABORT("DFT+U| An orbital magnetic quantum number was specified twice")
2400 : END IF
2401 : END IF
2402 : END DO
2403 : END DO
2404 : CALL section_vals_val_get(enforce_occupation_section, &
2405 : keyword_name="EPS_SCF", &
2406 4 : r_val=qs_kind%dft_plus_u%eps_scf)
2407 : CALL section_vals_val_get(enforce_occupation_section, &
2408 : keyword_name="MAX_SCF", &
2409 4 : i_val=i)
2410 4 : qs_kind%dft_plus_u%max_scf = MAX(-1, i)
2411 : CALL section_vals_val_get(enforce_occupation_section, &
2412 : keyword_name="SMEAR", &
2413 4 : l_val=qs_kind%dft_plus_u%smear)
2414 : END IF ! subsection enabled
2415 :
2416 90 : NULLIFY (tensorial_projector_section)
2417 : tensorial_projector_section => section_vals_get_subs_vals(dft_plus_u_section, &
2418 90 : subsection_name="ATOMIC_PROJECTOR_TENSORIAL")
2419 : subsection_enabled = .FALSE.
2420 : CALL section_vals_val_get(tensorial_projector_section, &
2421 : keyword_name="_SECTION_PARAMETERS_", &
2422 90 : l_val=subsection_enabled)
2423 90 : IF (subsection_enabled) THEN
2424 : CALL section_vals_val_get(tensorial_projector_section, &
2425 : keyword_name="SHELL_CHARGE", &
2426 8 : r_vals=proj_shell_charge)
2427 40 : qs_kind%dft_plus_u%proj_shell_charge = proj_shell_charge
2428 : END IF
2429 :
2430 90 : NULLIFY (minimum_tracking_linear_response_section)
2431 : minimum_tracking_linear_response_section => section_vals_get_subs_vals(dft_plus_u_section, &
2432 90 : subsection_name="MINIMUM_TRACKING_LINEAR_RESPONSE")
2433 : subsection_enabled = .FALSE.
2434 : CALL section_vals_val_get(minimum_tracking_linear_response_section, &
2435 : keyword_name="_SECTION_PARAMETERS_", &
2436 90 : l_val=subsection_enabled)
2437 90 : IF (subsection_enabled) THEN
2438 16 : qs_kind%dft_plus_u%do_mtlr = .TRUE.
2439 : CALL section_vals_val_get(minimum_tracking_linear_response_section, &
2440 : keyword_name="INDEX_PERTURBED_ATOM", &
2441 16 : i_val=qs_kind%dft_plus_u%lr_atom)
2442 16 : IF (qs_kind%dft_plus_u%lr_atom < 1) THEN
2443 0 : CPABORT("MTLR requires an explicit positive INDEX_PERTURBED_ATOM.")
2444 : END IF
2445 16 : NULLIFY (tmp_perturbation_strength)
2446 : CALL section_vals_val_get(minimum_tracking_linear_response_section, &
2447 : keyword_name="PERTURBATION_STRENGTH", &
2448 16 : r_vals=tmp_perturbation_strength)
2449 16 : IF (.NOT. ASSOCIATED(tmp_perturbation_strength)) THEN
2450 0 : CPABORT("MTLR perturbation strengths could not be read.")
2451 : END IF
2452 16 : IF (SIZE(tmp_perturbation_strength) < 3 .OR. SIZE(tmp_perturbation_strength) > 9) THEN
2453 : CALL cp_abort(__LOCATION__, &
2454 0 : "PERTURBATION_STRENGTH must contain between 3 and 9 values.")
2455 : END IF
2456 176 : IF (MAXVAL(tmp_perturbation_strength) - MINVAL(tmp_perturbation_strength) <= &
2457 : 100.0_dp*EPSILON(1.0_dp)) THEN
2458 0 : CPABORT("MTLR perturbation strengths must contain distinct values.")
2459 : END IF
2460 16 : IF (ALLOCATED(qs_kind%dft_plus_u%perturbation_strength)) THEN
2461 0 : DEALLOCATE (qs_kind%dft_plus_u%perturbation_strength)
2462 : END IF
2463 48 : ALLOCATE (qs_kind%dft_plus_u%perturbation_strength(SIZE(tmp_perturbation_strength)))
2464 96 : qs_kind%dft_plus_u%perturbation_strength(:) = tmp_perturbation_strength(:)
2465 : ! section_vals_val_get returns a pointer to storage owned by the
2466 : ! input tree. Only the copied component above is owned here.
2467 16 : NULLIFY (tmp_perturbation_strength)
2468 : END IF
2469 : END IF ! section enabled
2470 :
2471 : END IF
2472 :
2473 : ! Allocate and initialise the orbital basis set data set structure
2474 17311 : CALL init_orbital_pointers(5) ! debug the SUN optimizer
2475 :
2476 : ! BASIS and POTENTIAL read only when strictly necessary otherwise, even if not used
2477 : ! we just print misleading informations
2478 17311 : explicit_basis = .FALSE.
2479 17311 : IF (k_rep > 0) THEN
2480 : basis_section => section_vals_get_subs_vals(kind_section, "BASIS", i_rep_section=k_rep, &
2481 12045 : can_return_null=.TRUE.)
2482 12045 : CALL section_vals_get(basis_section, explicit=explicit_basis)
2483 : END IF
2484 :
2485 17311 : explicit_potential = .FALSE.
2486 17311 : IF (k_rep > 0) THEN
2487 : potential_section => section_vals_get_subs_vals(kind_section, "POTENTIAL", &
2488 12045 : i_rep_section=k_rep, can_return_null=.TRUE.)
2489 12045 : CALL section_vals_get(potential_section, explicit=explicit_potential)
2490 : END IF
2491 :
2492 17311 : explicit_kgpot = .FALSE.
2493 17311 : IF (k_rep > 0) THEN
2494 : kgpot_section => section_vals_get_subs_vals(kind_section, "KG_POTENTIAL", &
2495 12045 : i_rep_section=k_rep, can_return_null=.TRUE.)
2496 12045 : CALL section_vals_get(kgpot_section, explicit=explicit_kgpot)
2497 : END IF
2498 :
2499 17311 : CALL set_potential_file_names(potential_file_names, dft_potential_file_names, potential_fn_kind)
2500 :
2501 19547 : SELECT CASE (method_id)
2502 : CASE (do_method_rm1, do_method_am1, do_method_mndo, do_method_pdg, do_method_pm3, do_method_pm6, &
2503 : do_method_pm6fm, do_method_mndod, do_method_pnnl)
2504 : ! Allocate all_potential
2505 2236 : CALL allocate_potential(qs_kind%all_potential)
2506 2236 : CALL set_default_all_potential(qs_kind%all_potential, z, zeff_correction)
2507 2236 : CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
2508 2236 : IF (.NOT. ASSOCIATED(elec_conf)) THEN
2509 2236 : CALL get_potential(potential=qs_kind%all_potential, elec_conf=elec_conf)
2510 2236 : CALL set_qs_kind(qs_kind, elec_conf=elec_conf)
2511 : END IF
2512 2236 : CPASSERT(.NOT. qs_kind%floating)
2513 2236 : IF (qs_kind%ghost) THEN
2514 0 : CALL get_qs_kind(qs_kind=qs_kind, elec_conf=elec_conf)
2515 0 : elec_conf(:) = 0
2516 : CALL get_potential(potential=qs_kind%all_potential, &
2517 0 : elec_conf=elec_conf)
2518 0 : elec_conf(:) = 0
2519 : CALL set_potential(potential=qs_kind%all_potential, &
2520 : zeff=0.0_dp, &
2521 0 : zeff_correction=0.0_dp)
2522 : END IF
2523 :
2524 : ! Basis set (Parameters)
2525 : ! Setup proper semiempirical parameters
2526 2236 : check = .NOT. ASSOCIATED(qs_kind%se_parameter)
2527 2236 : CPASSERT(check)
2528 2236 : CALL semi_empirical_create(qs_kind%se_parameter)
2529 : ! Check if we allow p-orbitals on H
2530 438 : SELECT CASE (z)
2531 : CASE (1)
2532 2236 : IF (k_rep > 0) THEN
2533 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
2534 52 : keyword_name="SE_P_ORBITALS_ON_H", l_val=qs_kind%se_parameter%p_orbitals_on_h)
2535 : END IF
2536 : CASE DEFAULT
2537 : ! No special cases for other elements..
2538 : END SELECT
2539 : ! Set default parameters
2540 2236 : CALL section_vals_val_get(dft_section, "QS%SE%STO_NG", i_val=ngauss)
2541 2236 : CALL se_param_set_default(qs_kind%se_parameter, z, method_id)
2542 2236 : NULLIFY (tmp_basis_set)
2543 2236 : CALL init_se_param(qs_kind%se_parameter, tmp_basis_set, ngauss)
2544 2236 : CALL add_basis_set_to_container(qs_kind%basis_sets, tmp_basis_set, "ORB")
2545 : CALL init_potential(qs_kind%all_potential, itype="BARE", &
2546 2236 : zeff=qs_kind%se_parameter%zeff, zeff_correction=zeff_correction)
2547 2236 : qs_kind%se_parameter%zeff = qs_kind%se_parameter%zeff - zeff_correction
2548 :
2549 2236 : check = ((potential_name /= '') .OR. explicit_potential) .AND. .NOT. silent
2550 : IF (check) THEN
2551 : CALL cp_warn(__LOCATION__, &
2552 : "Information provided in the input file regarding POTENTIAL for KIND <"// &
2553 80 : TRIM(qs_kind%name)//"> will be ignored!")
2554 : END IF
2555 :
2556 2236 : check = ((k_rep > 0) .OR. explicit_basis) .AND. .NOT. silent
2557 630 : IF (check) THEN
2558 : CALL cp_warn(__LOCATION__, &
2559 : "Information provided in the input file regarding BASIS for KIND <"// &
2560 116 : TRIM(qs_kind%name)//"> will be ignored!")
2561 : END IF
2562 :
2563 : CASE (do_method_dftb)
2564 : ! Allocate all_potential
2565 630 : CALL allocate_potential(qs_kind%all_potential)
2566 630 : CALL set_default_all_potential(qs_kind%all_potential, z, zeff_correction)
2567 630 : CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
2568 630 : IF (.NOT. ASSOCIATED(elec_conf)) THEN
2569 630 : CALL get_potential(potential=qs_kind%all_potential, elec_conf=elec_conf)
2570 630 : CALL set_qs_kind(qs_kind, elec_conf=elec_conf)
2571 : END IF
2572 630 : CPASSERT(.NOT. qs_kind%floating)
2573 630 : IF (qs_kind%ghost) THEN
2574 0 : CALL get_qs_kind(qs_kind=qs_kind, elec_conf=elec_conf)
2575 0 : elec_conf(:) = 0
2576 : CALL get_potential(potential=qs_kind%all_potential, &
2577 0 : elec_conf=elec_conf)
2578 0 : elec_conf(:) = 0
2579 : CALL set_potential(potential=qs_kind%all_potential, &
2580 : zeff=0.0_dp, &
2581 0 : zeff_correction=0.0_dp)
2582 : END IF
2583 :
2584 630 : check = ((potential_name /= '') .OR. explicit_potential) .AND. .NOT. silent
2585 : IF (check) THEN
2586 : CALL cp_warn(__LOCATION__, &
2587 : "Information provided in the input file regarding POTENTIAL for KIND <"// &
2588 0 : TRIM(qs_kind%name)//"> will be ignored!")
2589 : END IF
2590 :
2591 630 : check = ((k_rep > 0) .OR. explicit_basis) .AND. .NOT. silent
2592 2614 : IF (check) THEN
2593 : CALL cp_warn(__LOCATION__, &
2594 : "Information provided in the input file regarding BASIS for KIND <"// &
2595 90 : TRIM(qs_kind%name)//"> will be ignored!")
2596 : END IF
2597 :
2598 : CASE (do_method_xtb)
2599 : ! Allocate all_potential
2600 2614 : CALL allocate_potential(qs_kind%all_potential)
2601 2614 : CALL set_default_all_potential(qs_kind%all_potential, z, zeff_correction)
2602 2614 : CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
2603 2614 : IF (.NOT. ASSOCIATED(elec_conf)) THEN
2604 2614 : CALL get_potential(potential=qs_kind%all_potential, elec_conf=elec_conf)
2605 2614 : CALL set_qs_kind(qs_kind, elec_conf=elec_conf)
2606 : END IF
2607 2614 : CPASSERT(.NOT. qs_kind%floating)
2608 2614 : IF (qs_kind%ghost) THEN
2609 0 : CALL get_qs_kind(qs_kind=qs_kind, elec_conf=elec_conf)
2610 0 : elec_conf(:) = 0
2611 : CALL get_potential(potential=qs_kind%all_potential, &
2612 0 : elec_conf=elec_conf)
2613 0 : elec_conf(:) = 0
2614 : CALL set_potential(potential=qs_kind%all_potential, &
2615 : zeff=0.0_dp, &
2616 0 : zeff_correction=0.0_dp)
2617 : END IF
2618 :
2619 2614 : check = ((potential_name /= '') .OR. explicit_potential) .AND. .NOT. silent
2620 : IF (check) THEN
2621 : CALL cp_warn(__LOCATION__, &
2622 : "Information provided in the input file regarding POTENTIAL for KIND <"// &
2623 0 : TRIM(qs_kind%name)//"> will be ignored!")
2624 : END IF
2625 :
2626 2614 : check = ((k_rep > 0) .OR. explicit_basis) .AND. .NOT. silent
2627 30 : IF (check) THEN
2628 : CALL cp_warn(__LOCATION__, &
2629 : "Information provided in the input file regarding BASIS for KIND <"// &
2630 0 : TRIM(qs_kind%name)//"> will be ignored!")
2631 : END IF
2632 :
2633 : CASE (do_method_pw)
2634 : ! PW DFT
2635 : ! Allocate and initialise the potential data set structure
2636 30 : IF (potential_name /= '') THEN
2637 30 : SELECT CASE (TRIM(potential_type))
2638 : CASE ("ALL", "ECP", "CNEO")
2639 : CALL cp_abort(__LOCATION__, &
2640 : "PW DFT calculations only with potential type UPF or GTH possible."// &
2641 : " <"//TRIM(potential_type)//"> was specified "// &
2642 0 : "for the atomic kind <"//TRIM(qs_kind%name))
2643 : CASE ("GTH")
2644 8 : CALL allocate_potential(qs_kind%gth_potential)
2645 : CALL read_potential(qs_kind%element_symbol, potential_name, &
2646 : qs_kind%gth_potential, zeff_correction, para_env, &
2647 : potential_file_names, potential_section, update_input, &
2648 8 : monovalent=qs_kind%monovalent)
2649 8 : CALL set_potential(qs_kind%gth_potential, z=z)
2650 8 : CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
2651 8 : IF (.NOT. ASSOCIATED(elec_conf)) THEN
2652 8 : CALL get_potential(potential=qs_kind%gth_potential, elec_conf=elec_conf)
2653 8 : CALL set_qs_kind(qs_kind, elec_conf=elec_conf)
2654 : ELSE
2655 0 : CALL set_potential(potential=qs_kind%gth_potential, elec_conf=elec_conf)
2656 : END IF
2657 : CASE ("UPF")
2658 2310 : ALLOCATE (qs_kind%upf_potential)
2659 22 : qs_kind%upf_potential%zion = 0
2660 22 : qs_kind%upf_potential%filename = ADJUSTL(TRIM(potential_name))
2661 22 : CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
2662 22 : IF (.NOT. ASSOCIATED(elec_conf)) THEN
2663 22 : CALL set_qs_kind(qs_kind, elec_conf=qs_kind%upf_potential%econf)
2664 : END IF
2665 : CASE DEFAULT
2666 : CALL cp_abort(__LOCATION__, &
2667 : "An invalid potential type <"// &
2668 : TRIM(potential_type)//"> was specified "// &
2669 : "for the atomic kind <"// &
2670 30 : TRIM(qs_kind%name))
2671 : END SELECT
2672 : ELSE
2673 : CALL cp_abort(__LOCATION__, &
2674 : "No potential type was defined for the "// &
2675 0 : "atomic kind <"//TRIM(qs_kind%name)//">")
2676 : END IF
2677 :
2678 : CASE DEFAULT
2679 :
2680 : ! set ngauss for STO expansion
2681 11801 : CALL section_vals_val_get(dft_section, "QS%STO_NG", i_val=ngauss)
2682 : ! Allocate and initialise the basis set data set structure
2683 : ! first external basis sets
2684 25959 : DO i = 1, nb_rep
2685 28312 : SELECT CASE (basis_set_form(i))
2686 : CASE ("GTO")
2687 14154 : NULLIFY (tmp_basis_set)
2688 14154 : CALL allocate_gto_basis_set(tmp_basis_set)
2689 : CALL read_gto_basis_set(qs_kind%element_symbol, basis_set_name(i), &
2690 14154 : tmp_basis_set, para_env, dft_section)
2691 : CASE ("STO")
2692 4 : NULLIFY (sto_basis_set)
2693 4 : CALL allocate_sto_basis_set(sto_basis_set)
2694 : CALL read_sto_basis_set(qs_kind%element_symbol, basis_set_name(i), &
2695 4 : sto_basis_set, para_env, dft_section)
2696 4 : NULLIFY (tmp_basis_set)
2697 4 : CALL create_gto_from_sto_basis(sto_basis_set, tmp_basis_set, ngauss)
2698 4 : CALL deallocate_sto_basis_set(sto_basis_set)
2699 : CASE DEFAULT
2700 : CALL cp_abort(__LOCATION__, &
2701 : "Invalid basis set form "//TRIM(basis_set_form(i))// &
2702 14158 : "for atomic kind <"//TRIM(qs_kind%name)//">")
2703 : END SELECT
2704 14158 : tmp = basis_set_type(i)
2705 14158 : CALL uppercase(tmp)
2706 25959 : CALL add_basis_set_to_container(qs_kind%basis_sets, tmp_basis_set, tmp)
2707 : END DO
2708 : ! now explicit basis sets
2709 11801 : IF (explicit_basis) THEN
2710 158 : CALL section_vals_get(basis_section, n_repetition=nexp)
2711 316 : DO i = 1, nexp
2712 158 : NULLIFY (tmp_basis_set)
2713 158 : CALL allocate_gto_basis_set(tmp_basis_set)
2714 : CALL read_gto_basis_set(qs_kind%element_symbol, basis_type, &
2715 158 : tmp_basis_set, basis_section, i, dft_section)
2716 158 : tmp = basis_type
2717 158 : CALL uppercase(tmp)
2718 316 : CALL add_basis_set_to_container(qs_kind%basis_sets, tmp_basis_set, tmp)
2719 : END DO
2720 : END IF
2721 : ! combine multiple basis sets
2722 247821 : DO i = 1, SIZE(qs_kind%basis_sets)
2723 236020 : NULLIFY (tmp_basis_set)
2724 : CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis_set, &
2725 236020 : inumbas=i, basis_type=basis_type)
2726 236020 : IF (basis_type == "") CYCLE
2727 14316 : jj = i
2728 283675 : DO j = i + 1, SIZE(qs_kind%basis_sets)
2729 269359 : jj = jj + 1
2730 269359 : NULLIFY (sup_basis_set)
2731 : CALL get_basis_from_container(qs_kind%basis_sets, basis_set=sup_basis_set, &
2732 269359 : inumbas=jj, basis_type=tmp)
2733 283675 : IF (basis_type == tmp) THEN
2734 : ! we found a match, combine the basis sets and delete the second
2735 0 : CALL combine_basis_sets(tmp_basis_set, sup_basis_set)
2736 0 : CALL remove_basis_from_container(qs_kind%basis_sets, jj)
2737 0 : jj = jj - 1
2738 : END IF
2739 : END DO
2740 247821 : NULLIFY (sup_basis_set)
2741 : END DO
2742 :
2743 : ! check that we have an orbital basis set
2744 11801 : nobasis = .TRUE.
2745 247821 : DO i = 1, SIZE(qs_kind%basis_sets)
2746 236020 : NULLIFY (tmp_basis_set)
2747 : CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis_set, &
2748 236020 : inumbas=i, basis_type=basis_type)
2749 247821 : IF (basis_type == "ORB") nobasis = .FALSE.
2750 : END DO
2751 11801 : IF (nobasis) THEN
2752 : CALL cp_abort(__LOCATION__, &
2753 : "No basis set type was defined for the "// &
2754 0 : "atomic kind <"//TRIM(qs_kind%name)//">")
2755 : END IF
2756 :
2757 : ! If Ghost atom we don't need to allocate/initialize anything connected to POTENTIAL
2758 11801 : IF (qs_kind%ghost .OR. qs_kind%floating) THEN
2759 152 : IF (ASSOCIATED(qs_kind%elec_conf)) qs_kind%elec_conf = 0
2760 : ELSE
2761 : ! Allocate and initialise the potential data set structure
2762 11649 : IF ((potential_name /= '') .OR. explicit_potential) THEN
2763 13119 : SELECT CASE (TRIM(potential_type))
2764 : CASE ("ALL")
2765 1470 : CALL allocate_potential(qs_kind%all_potential)
2766 : CALL read_potential(qs_kind%element_symbol, potential_name, &
2767 : qs_kind%all_potential, zeff_correction, para_env, &
2768 1470 : potential_file_names, potential_section, update_input)
2769 1470 : CALL set_potential(qs_kind%all_potential, z=z)
2770 1470 : CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
2771 1470 : IF (.NOT. ASSOCIATED(elec_conf)) THEN
2772 1470 : CALL get_potential(potential=qs_kind%all_potential, elec_conf=elec_conf)
2773 1470 : CALL set_qs_kind(qs_kind, elec_conf=elec_conf)
2774 : ELSE
2775 0 : CALL set_potential(potential=qs_kind%all_potential, elec_conf=elec_conf)
2776 : END IF
2777 : CASE ("GTH")
2778 10065 : CALL allocate_potential(qs_kind%gth_potential)
2779 : CALL read_potential(qs_kind%element_symbol, potential_name, &
2780 : qs_kind%gth_potential, zeff_correction, para_env, &
2781 : potential_file_names, potential_section, update_input, &
2782 10065 : monovalent=qs_kind%monovalent)
2783 10065 : CALL set_potential(qs_kind%gth_potential, z=z)
2784 10065 : CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
2785 10065 : IF (.NOT. ASSOCIATED(elec_conf)) THEN
2786 10065 : CALL get_potential(potential=qs_kind%gth_potential, elec_conf=elec_conf)
2787 10065 : CALL set_qs_kind(qs_kind, elec_conf=elec_conf)
2788 : ELSE
2789 0 : CALL set_potential(potential=qs_kind%gth_potential, elec_conf=elec_conf)
2790 : END IF
2791 : CASE ("ECP")
2792 82 : CALL allocate_potential(qs_kind%sgp_potential)
2793 82 : CALL get_potential(qs_kind%sgp_potential, description=description)
2794 : CALL read_ecp_potential(ptable(z)%symbol, ecppot, &
2795 82 : potential_name, potential_file_names, potential_section)
2796 82 : IF (ecp_semi_local) THEN
2797 82 : description(1) = "Semi-local Gaussian pseudopotential "
2798 82 : description(2) = "ECP "//TRIM(potential_name)
2799 82 : description(3) = "LIBGRPP: A. V. Oleynichenko et al., Symmetry 15 197 2023"
2800 : description(4) = " "
2801 : ELSE
2802 0 : description(4) = "ECP "//TRIM(potential_name)
2803 : END IF
2804 : CALL set_potential(qs_kind%sgp_potential, name=ecppot%pname, description=description, &
2805 : zeff=ecppot%zion, z=z, ecp_local=.TRUE., ecp_semi_local=ecp_semi_local, &
2806 : nloc=ecppot%nloc, nrloc=ecppot%nrloc, aloc=ecppot%aloc, bloc=ecppot%bloc, &
2807 82 : has_nlcc=.FALSE.)
2808 : CALL set_potential(qs_kind%sgp_potential, sl_lmax=ecppot%lmax, &
2809 82 : npot=ecppot%npot, nrpot=ecppot%nrpot, apot=ecppot%apot, bpot=ecppot%bpot)
2810 : ! convert PP
2811 82 : IF (.NOT. ecp_semi_local) THEN
2812 0 : CPABORT("ECPs are only well tested in their semi-local form")
2813 0 : CALL get_qs_kind(qs_kind, basis_set=orb_basis_set)
2814 0 : CALL sgp_construction(sgp_pot=sgppot, ecp_pot=ecppot, orb_basis=orb_basis_set, error=error)
2815 0 : IF (iounit > 0 .AND. .NOT. silent) THEN
2816 0 : WRITE (iounit, "(/,T2,'PP Transformation for ',A)") TRIM(ecppot%pname)
2817 0 : IF (sgppot%has_local) THEN
2818 0 : WRITE (iounit, "(T8,'Accuracy for local part:',T41,F10.3,'%',T61,F20.12)") error(4), error(1)
2819 : END IF
2820 0 : IF (sgppot%has_nonlocal) THEN
2821 0 : WRITE (iounit, "(T8,'Accuracy for nonlocal part:',T41,F10.3,'%',T61,F20.12)") error(5), error(2)
2822 : END IF
2823 0 : IF (sgppot%has_nlcc) THEN
2824 0 : WRITE (iounit, "(T8,'Accuracy for NLCC density:',T61,F20.12)") error(3)
2825 : END IF
2826 : END IF
2827 : END IF
2828 82 : IF (sgppot%has_nonlocal) THEN
2829 : CALL set_potential(qs_kind%sgp_potential, n_nonlocal=sgppot%n_nonlocal, lmax=sgppot%lmax, &
2830 0 : is_nonlocal=sgppot%is_nonlocal)
2831 0 : nnl = sgppot%n_nonlocal
2832 0 : nppnl = 0
2833 0 : DO l = 0, sgppot%lmax
2834 0 : nppnl = nppnl + nnl*nco(l)
2835 : END DO
2836 0 : l = sgppot%lmax
2837 0 : ALLOCATE (a_nl(nnl), h_nl(nnl, 0:l), c_nl(nnl, nnl, 0:l))
2838 0 : a_nl(:) = sgppot%a_nonlocal(:)
2839 0 : h_nl(:, :) = sgppot%h_nonlocal(:, :)
2840 0 : c_nl(:, :, :) = sgppot%c_nonlocal(:, :, :)
2841 0 : CALL set_potential(qs_kind%sgp_potential, nppnl=nppnl, a_nonlocal=a_nl, h_nonlocal=h_nl, c_nonlocal=c_nl)
2842 : ELSE
2843 82 : CALL set_potential(qs_kind%sgp_potential, n_nonlocal=0, lmax=-1, is_nonlocal=sgppot%is_nonlocal)
2844 82 : CALL set_potential(qs_kind%sgp_potential, nppnl=0)
2845 : END IF
2846 : !
2847 82 : CPASSERT(.NOT. sgppot%has_local)
2848 82 : CPASSERT(.NOT. sgppot%has_nlcc)
2849 : ! core
2850 82 : rc = 0.5_dp*qs_kind%covalent_radius*angstrom
2851 82 : rc = MAX(rc, 0.2_dp)
2852 82 : rc = MIN(rc, 1.0_dp)
2853 82 : alpha = 1.0_dp/(2.0_dp*rc**2)
2854 82 : ccore = ecppot%zion*SQRT((alpha/pi)**3)
2855 : CALL set_potential(qs_kind%sgp_potential, alpha_core_charge=alpha, ccore_charge=ccore, &
2856 82 : core_charge_radius=rc)
2857 82 : CALL atom_sgp_release(sgppot)
2858 82 : CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
2859 82 : IF (.NOT. ASSOCIATED(elec_conf)) THEN
2860 82 : CALL set_qs_kind(qs_kind, elec_conf=ecppot%econf)
2861 : END IF
2862 82 : CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
2863 82 : CALL set_potential(qs_kind%sgp_potential, elec_conf=elec_conf)
2864 : CASE ("UPF")
2865 24 : CALL allocate_potential(qs_kind%sgp_potential)
2866 24 : CALL get_potential(qs_kind%sgp_potential, description=description)
2867 24 : description(4) = "UPF "//TRIM(potential_name)
2868 24 : CALL atom_read_upf(upfpot, potential_name)
2869 24 : IF (TRIM(upfpot%pseudo_type) == "NC") THEN
2870 24 : ALLOCATE (qs_kind%sgp_potential%upf_data, SOURCE=upfpot)
2871 : END IF
2872 : CALL set_potential(qs_kind%sgp_potential, name=upfpot%pname, description=description, &
2873 24 : zeff=upfpot%zion, z=z, has_nlcc=upfpot%core_correction)
2874 : ! convert pp
2875 24 : CALL sgp_construction(sgp_pot=sgppot, upf_pot=upfpot, error=error)
2876 24 : IF (iounit > 0 .AND. .NOT. silent) THEN
2877 12 : WRITE (iounit, "(/,T2,'PP Transformation for ',A)") TRIM(upfpot%pname)
2878 12 : IF (sgppot%has_local) THEN
2879 12 : WRITE (iounit, "(T8,'Auxiliary local Gaussian fit error:',T61,F20.12)") error(1)
2880 : END IF
2881 12 : IF (sgppot%has_nonlocal) THEN
2882 6 : IF (ALLOCATED(qs_kind%sgp_potential%upf_data)) THEN
2883 6 : WRITE (iounit, "(T8,'Unused nonlocal Gaussian fit error:',T61,F20.12)") error(2)
2884 : ELSE
2885 0 : WRITE (iounit, "(T8,'Accuracy for nonlocal part:',T61,F20.12)") error(2)
2886 : END IF
2887 : END IF
2888 12 : IF (sgppot%has_nlcc) THEN
2889 3 : WRITE (iounit, "(T8,'Accuracy for NLCC density:',T61,F20.12)") error(3)
2890 : END IF
2891 : END IF
2892 24 : IF (sgppot%has_nonlocal) THEN
2893 : CALL set_potential(qs_kind%sgp_potential, n_nonlocal=sgppot%n_nonlocal, lmax=sgppot%lmax, &
2894 12 : is_nonlocal=sgppot%is_nonlocal)
2895 12 : nnl = sgppot%n_nonlocal
2896 12 : nppnl = 0
2897 36 : DO l = 0, sgppot%lmax
2898 36 : nppnl = nppnl + nnl*nco(l)
2899 : END DO
2900 12 : l = sgppot%lmax
2901 120 : ALLOCATE (a_nl(nnl), h_nl(nnl, 0:l), c_nl(nnl, nnl, 0:l))
2902 108 : a_nl(:) = sgppot%a_nonlocal(:)
2903 228 : h_nl(:, :) = sgppot%h_nonlocal(:, :)
2904 1764 : c_nl(:, :, :) = sgppot%c_nonlocal(:, :, :)
2905 12 : CALL set_potential(qs_kind%sgp_potential, nppnl=nppnl, a_nonlocal=a_nl, h_nonlocal=h_nl, c_nonlocal=c_nl)
2906 : ELSE
2907 12 : CALL set_potential(qs_kind%sgp_potential, n_nonlocal=0, lmax=-1, is_nonlocal=sgppot%is_nonlocal)
2908 12 : CALL set_potential(qs_kind%sgp_potential, nppnl=0)
2909 : END IF
2910 24 : CPASSERT(sgppot%has_local)
2911 : ! core
2912 24 : rc = sgppot%ac_local
2913 24 : alpha = 1.0_dp/(2.0_dp*rc**2)
2914 24 : ccore = upfpot%zion*SQRT((alpha/pi)**3)
2915 : CALL set_potential(qs_kind%sgp_potential, alpha_core_charge=alpha, ccore_charge=ccore, &
2916 24 : core_charge_radius=rc)
2917 : ! local potential
2918 24 : nloc = sgppot%n_local
2919 96 : ALLOCATE (aloc(nloc), cloc(nloc))
2920 312 : aloc(1:nloc) = sgppot%a_local(1:nloc)
2921 312 : cloc(1:nloc) = sgppot%c_local(1:nloc)
2922 24 : CALL set_potential(qs_kind%sgp_potential, n_local=nloc, a_local=aloc, c_local=cloc)
2923 24 : IF (sgppot%has_nlcc) THEN
2924 6 : nlcc = sgppot%n_nlcc
2925 24 : ALLOCATE (anlcc(nlcc), cnlcc(nlcc))
2926 294 : anlcc(1:nlcc) = sgppot%a_nlcc(1:nlcc)
2927 294 : cnlcc(1:nlcc) = sgppot%c_nlcc(1:nlcc)
2928 6 : CALL set_potential(qs_kind%sgp_potential, has_nlcc=.TRUE., n_nlcc=nlcc, a_nlcc=anlcc, c_nlcc=cnlcc)
2929 : END IF
2930 24 : CALL set_potential(qs_kind%sgp_potential, z=z)
2931 24 : CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
2932 24 : IF (.NOT. ASSOCIATED(elec_conf)) THEN
2933 24 : CALL set_qs_kind(qs_kind, elec_conf=upfpot%econf)
2934 : END IF
2935 24 : CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
2936 24 : CALL set_potential(qs_kind%sgp_potential, elec_conf=elec_conf)
2937 24 : CALL atom_release_upf(upfpot)
2938 24 : CALL atom_sgp_release(sgppot)
2939 : CASE ("CNEO")
2940 8 : IF (zeff_correction /= 0.0_dp) THEN
2941 0 : CPABORT("CORE_CORRECTION is not compatible with CNEO")
2942 : END IF
2943 8 : CALL allocate_cneo_potential(qs_kind%cneo_potential)
2944 8 : CALL set_cneo_potential(qs_kind%cneo_potential, z=z)
2945 8 : mass = 0.0_dp
2946 : ! Input mass is the mass of the neutral atom, not the nucleus.
2947 : ! The mass of electrons will get subtracted later.
2948 : ! In principle, the most abundant pure isotope mass should be used.
2949 8 : IF (k_rep > 0) THEN
2950 : CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
2951 8 : keyword_name="MASS", n_rep_val=i)
2952 8 : IF (i > 0) CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
2953 2 : keyword_name="MASS", r_val=mass)
2954 : END IF
2955 : ! Remove electron mass from atomic mass to get nuclear mass
2956 8 : IF (mass - REAL(z, dp)*0.000548579909_dp > 0.0_dp) THEN
2957 2 : mass = mass - REAL(z, dp)*0.000548579909_dp
2958 2 : CALL set_cneo_potential(qs_kind%cneo_potential, mass=mass)
2959 : END IF
2960 : ! In case the mass is not set by user, get the default mass from z
2961 8 : CALL get_cneo_potential(qs_kind%cneo_potential, mass=mass)
2962 : ! Warn if the mass is from ptable
2963 8 : IF (ABS(mass + REAL(z, dp)*0.000548579909_dp - ptable(z)%amass) < 1.e-4_dp) THEN
2964 : CALL cp_warn(__LOCATION__, &
2965 : "Atomic mass of the atomic kind <"//TRIM(qs_kind%name)// &
2966 : "> is very close to its average mass. Is it a pure isotope? "// &
2967 : "Pure isotopes are preferable for CNEO. "// &
2968 0 : "(e.g., mass of 1H is 1.007825, not 1.00794)")
2969 : END IF
2970 8 : CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
2971 8 : IF (.NOT. ASSOCIATED(elec_conf)) THEN
2972 8 : CALL get_cneo_potential(potential=qs_kind%cneo_potential, elec_conf=elec_conf)
2973 8 : CALL set_qs_kind(qs_kind, elec_conf=elec_conf)
2974 : ELSE
2975 0 : CALL set_cneo_potential(potential=qs_kind%cneo_potential, elec_conf=elec_conf)
2976 : END IF
2977 : CASE DEFAULT
2978 : CALL cp_abort(__LOCATION__, &
2979 : "An invalid potential type <"// &
2980 : TRIM(potential_name)//"> was specified "// &
2981 : "for the atomic kind <"// &
2982 11657 : TRIM(qs_kind%name))
2983 : END SELECT
2984 : ELSE
2985 : CALL cp_abort(__LOCATION__, &
2986 : "No potential type was defined for the "// &
2987 0 : "atomic kind <"//TRIM(qs_kind%name)//">")
2988 : END IF
2989 :
2990 11649 : CALL check_potential_basis_compatibility(qs_kind)
2991 :
2992 : ! Allocate and initialise the potential data set structure
2993 11649 : IF ((kgpot_name /= '') .OR. explicit_kgpot) THEN
2994 11649 : ipos = INDEX(kgpot_name, "-")
2995 11649 : IF (ipos > 1) THEN
2996 24 : kgpot_type = kgpot_name(:ipos - 1)
2997 : ELSE
2998 11625 : kgpot_type = kgpot_name
2999 : END IF
3000 11649 : CALL uppercase(kgpot_type)
3001 :
3002 11673 : SELECT CASE (TRIM(kgpot_type))
3003 : CASE ("TNADD")
3004 : CALL set_potential_file_names(potential_file_names, dft_potential_file_names, &
3005 24 : kg_potential_fn_kind)
3006 24 : CALL allocate_potential(qs_kind%tnadd_potential)
3007 : CALL read_potential(qs_kind%element_symbol, kgpot_name, &
3008 : qs_kind%tnadd_potential, para_env, &
3009 24 : potential_file_names, kgpot_section, update_input)
3010 : CASE ("NONE")
3011 11625 : NULLIFY (qs_kind%tnadd_potential)
3012 : CASE DEFAULT
3013 : CALL cp_abort(__LOCATION__, &
3014 : "An invalid kg_potential type <"// &
3015 : TRIM(potential_name)//"> was specified "// &
3016 : "for the atomic kind <"// &
3017 11649 : TRIM(qs_kind%name))
3018 : END SELECT
3019 : END IF
3020 : END IF
3021 :
3022 : ! check that we have a nuclear orbital basis set if CNEO is requested
3023 11801 : nobasis_nuc = ASSOCIATED(qs_kind%cneo_potential)
3024 247821 : DO i = 1, SIZE(qs_kind%basis_sets)
3025 236020 : NULLIFY (tmp_basis_set)
3026 : CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis_set, &
3027 236020 : inumbas=i, basis_type=basis_type)
3028 247821 : IF (basis_type == "NUC") THEN
3029 8 : nobasis_nuc = .FALSE.
3030 8 : IF (.NOT. ASSOCIATED(qs_kind%cneo_potential)) THEN
3031 : CALL cp_warn(__LOCATION__, &
3032 : "POTENTIAL is not set to CNEO, NUC type basis set for KIND <"// &
3033 0 : TRIM(qs_kind%name)//"> will be ignored!")
3034 : END IF
3035 : END IF
3036 : END DO
3037 31348 : IF (nobasis_nuc) THEN
3038 : CALL cp_abort(__LOCATION__, &
3039 : "No NUC type basis set was defined for the "// &
3040 : "atomic kind <"//TRIM(qs_kind%name)//">, which is required by "// &
3041 0 : "POTENTIAL CNEO.")
3042 : END IF
3043 : END SELECT
3044 :
3045 17311 : DEALLOCATE (dft_potential_file_names, potential_file_names)
3046 17311 : CALL timestop(handle)
3047 :
3048 10421222 : END SUBROUTINE read_qs_kind
3049 :
3050 : ! **************************************************************************************************
3051 : !> \brief Ensure pseudo-potential and basis set were optimized for same number of valence electrons
3052 : !> \param qs_kind ...
3053 : !> \author Ole Schuett
3054 : ! **************************************************************************************************
3055 11649 : SUBROUTINE check_potential_basis_compatibility(qs_kind)
3056 : TYPE(qs_kind_type), INTENT(INOUT) :: qs_kind
3057 :
3058 : CHARACTER(LEN=default_string_length) :: name
3059 : INTEGER :: nbs, npp
3060 : TYPE(gth_potential_type), POINTER :: gth_potential
3061 : TYPE(gto_basis_set_type), POINTER :: basis_set
3062 :
3063 11649 : CALL get_qs_kind(qs_kind, name=name, gth_potential=gth_potential, basis_set=basis_set)
3064 :
3065 11649 : npp = -1; nbs = -1
3066 11649 : IF (ASSOCIATED(gth_potential)) THEN
3067 10065 : npp = parse_valence_electrons(gth_potential%aliases)
3068 : END IF
3069 11649 : IF (ASSOCIATED(basis_set)) THEN
3070 11649 : nbs = parse_valence_electrons(basis_set%aliases)
3071 : END IF
3072 :
3073 11649 : IF (npp >= 0 .AND. nbs >= 0 .AND. npp /= nbs) THEN
3074 : CALL cp_abort(__LOCATION__, "Basis-set and pseudo-potential of atomic kind '"//TRIM(name)//"'"// &
3075 0 : " were optimized for different valence electron numbers.")
3076 : END IF
3077 :
3078 11649 : END SUBROUTINE check_potential_basis_compatibility
3079 :
3080 : ! **************************************************************************************************
3081 : !> \brief Tries to parse valence eletron number using "-QXXX" notation, returns -1 if not found.
3082 : !> \param string ...
3083 : !> \return ...
3084 : !> \author Ole Schuett
3085 : ! **************************************************************************************************
3086 21714 : FUNCTION parse_valence_electrons(string) RESULT(n)
3087 : CHARACTER(*) :: string
3088 : INTEGER :: n
3089 :
3090 : INTEGER :: i, istat, j
3091 :
3092 21714 : i = INDEX(string, "-Q", .TRUE.)
3093 21714 : IF (i == 0) THEN
3094 7404 : n = -1
3095 : ELSE
3096 14310 : j = SCAN(string(i + 2:), "- ")
3097 14310 : READ (string(i + 2:i + j), '(I3)', iostat=istat) n
3098 14310 : IF (istat /= 0) n = -1
3099 : END IF
3100 :
3101 21714 : END FUNCTION parse_valence_electrons
3102 :
3103 : ! **************************************************************************************************
3104 : !> \brief Read an atomic kind set data set from the input file.
3105 : !> \param qs_kind_set ...
3106 : !> \param atomic_kind_set ...
3107 : !> \param kind_section ...
3108 : !> \param para_env ...
3109 : !> \param force_env_section ...
3110 : !> \param silent ...
3111 : ! **************************************************************************************************
3112 9376 : SUBROUTINE create_qs_kind_set(qs_kind_set, atomic_kind_set, kind_section, para_env, &
3113 : force_env_section, silent)
3114 :
3115 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3116 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3117 : TYPE(section_vals_type), POINTER :: kind_section
3118 : TYPE(mp_para_env_type), POINTER :: para_env
3119 : TYPE(section_vals_type), POINTER :: force_env_section
3120 : LOGICAL, INTENT(IN) :: silent
3121 :
3122 : CHARACTER(len=*), PARAMETER :: routineN = 'create_qs_kind_set'
3123 :
3124 : INTEGER :: handle, ikind, method, nkind, qs_method
3125 : LOGICAL :: no_fail
3126 :
3127 9376 : CALL timeset(routineN, handle)
3128 :
3129 9376 : IF (ASSOCIATED(qs_kind_set)) CPABORT("create_qs_kind_set: qs_kind_set already associated")
3130 9376 : IF (.NOT. ASSOCIATED(atomic_kind_set)) CPABORT("create_qs_kind_set: atomic_kind_set not associated")
3131 :
3132 9376 : no_fail = .FALSE.
3133 :
3134 : ! Between all methods only SE and DFTB/xTB may not need a KIND section.
3135 9376 : CALL section_vals_val_get(force_env_section, "METHOD", i_val=method)
3136 9376 : IF (method == do_qs) THEN
3137 9348 : CALL section_vals_val_get(force_env_section, "DFT%QS%METHOD", i_val=qs_method)
3138 998 : SELECT CASE (qs_method)
3139 : CASE (do_method_mndo, do_method_am1, do_method_pm3, do_method_pm6fm, do_method_pm6, &
3140 : do_method_pdg, do_method_rm1, do_method_mndod, do_method_pnnl)
3141 998 : no_fail = .TRUE.
3142 : CASE (do_method_dftb)
3143 298 : no_fail = .TRUE.
3144 : CASE (do_method_xtb)
3145 9348 : no_fail = .TRUE.
3146 : END SELECT
3147 28 : ELSE IF (method == do_sirius) THEN
3148 24 : qs_method = do_method_pw
3149 : ELSE
3150 4 : qs_method = method
3151 : END IF
3152 :
3153 9376 : nkind = SIZE(atomic_kind_set)
3154 232959 : ALLOCATE (qs_kind_set(nkind))
3155 :
3156 26687 : DO ikind = 1, nkind
3157 17311 : qs_kind_set(ikind)%name = atomic_kind_set(ikind)%name
3158 17311 : qs_kind_set(ikind)%element_symbol = atomic_kind_set(ikind)%element_symbol
3159 17311 : qs_kind_set(ikind)%natom = atomic_kind_set(ikind)%natom
3160 : CALL read_qs_kind(qs_kind_set(ikind), kind_section, para_env, force_env_section, &
3161 26687 : no_fail, qs_method, silent)
3162 : END DO
3163 :
3164 9376 : CALL timestop(handle)
3165 :
3166 18752 : END SUBROUTINE create_qs_kind_set
3167 :
3168 : ! **************************************************************************************************
3169 : !> \brief This routines should perform only checks. no settings are allowed at
3170 : !> this level anymore..
3171 : !> \param qs_kind ...
3172 : !> \param dft_control ...
3173 : !> \param subsys_section ...
3174 : ! **************************************************************************************************
3175 17225 : SUBROUTINE check_qs_kind(qs_kind, dft_control, subsys_section)
3176 :
3177 : TYPE(qs_kind_type), POINTER :: qs_kind
3178 : TYPE(dft_control_type), INTENT(IN) :: dft_control
3179 : TYPE(section_vals_type), POINTER :: subsys_section
3180 :
3181 : INTEGER :: gfn_type
3182 : LOGICAL :: defined
3183 : TYPE(qs_dftb_atom_type), POINTER :: dftb_parameter
3184 : TYPE(semi_empirical_type), POINTER :: se_parameter
3185 : TYPE(xtb_atom_type), POINTER :: xtb_parameter
3186 :
3187 17225 : IF (dft_control%qs_control%semi_empirical) THEN
3188 2236 : CALL get_qs_kind(qs_kind, se_parameter=se_parameter)
3189 2236 : CPASSERT(ASSOCIATED(se_parameter))
3190 2236 : CALL get_se_param(se_parameter, defined=defined)
3191 2236 : CPASSERT(defined)
3192 2236 : CALL write_se_param(se_parameter, subsys_section)
3193 14989 : ELSE IF (dft_control%qs_control%dftb) THEN
3194 630 : CALL get_qs_kind(qs_kind, dftb_parameter=dftb_parameter)
3195 630 : CPASSERT(ASSOCIATED(dftb_parameter))
3196 630 : CALL get_dftb_atom_param(dftb_parameter, defined=defined)
3197 630 : CPASSERT(defined)
3198 630 : CALL write_dftb_atom_param(dftb_parameter, subsys_section)
3199 14359 : ELSE IF (dft_control%qs_control%xtb) THEN
3200 2614 : IF (.NOT. (dft_control%qs_control%xtb_control%do_tblite)) THEN
3201 2280 : CALL get_qs_kind(qs_kind, xtb_parameter=xtb_parameter)
3202 2280 : CPASSERT(ASSOCIATED(xtb_parameter))
3203 2280 : gfn_type = dft_control%qs_control%xtb_control%gfn_type
3204 2280 : CALL write_xtb_atom_param(xtb_parameter, gfn_type, subsys_section)
3205 : END IF
3206 : END IF
3207 :
3208 17225 : END SUBROUTINE check_qs_kind
3209 :
3210 : ! **************************************************************************************************
3211 : !> \brief ...
3212 : !> \param qs_kind_set ...
3213 : !> \param dft_control ...
3214 : !> \param subsys_section ...
3215 : ! **************************************************************************************************
3216 9324 : SUBROUTINE check_qs_kind_set(qs_kind_set, dft_control, subsys_section)
3217 :
3218 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3219 : TYPE(dft_control_type), INTENT(IN) :: dft_control
3220 : TYPE(section_vals_type), POINTER :: subsys_section
3221 :
3222 : CHARACTER(len=*), PARAMETER :: routineN = 'check_qs_kind_set'
3223 :
3224 : INTEGER :: handle, ikind, nkind
3225 : TYPE(qs_kind_type), POINTER :: qs_kind
3226 :
3227 9324 : CALL timeset(routineN, handle)
3228 9324 : IF (ASSOCIATED(qs_kind_set)) THEN
3229 9324 : nkind = SIZE(qs_kind_set)
3230 26549 : DO ikind = 1, nkind
3231 17225 : qs_kind => qs_kind_set(ikind)
3232 26549 : CALL check_qs_kind(qs_kind, dft_control, subsys_section)
3233 : END DO
3234 9324 : IF (dft_control%qs_control%xtb) THEN
3235 : CALL write_xtb_kab_param(qs_kind_set, subsys_section, &
3236 1232 : dft_control%qs_control%xtb_control)
3237 : END IF
3238 : ELSE
3239 0 : CPABORT("The pointer qs_kind_set is not associated")
3240 : END IF
3241 9324 : CALL timestop(handle)
3242 9324 : END SUBROUTINE check_qs_kind_set
3243 :
3244 : ! **************************************************************************************************
3245 : !> \brief ...
3246 : !> \param qs_kind_set ...
3247 : !> \param subsys_section ...
3248 : !> \param xtb_control ...
3249 : ! **************************************************************************************************
3250 1232 : SUBROUTINE write_xtb_kab_param(qs_kind_set, subsys_section, xtb_control)
3251 :
3252 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3253 : TYPE(section_vals_type), POINTER :: subsys_section
3254 : TYPE(xtb_control_type), POINTER :: xtb_control
3255 :
3256 : CHARACTER(LEN=default_string_length) :: aname, bname
3257 : INTEGER :: ikind, io_unit, jkind, nkind, za, zb
3258 : TYPE(cp_logger_type), POINTER :: logger
3259 : TYPE(qs_kind_type), POINTER :: qs_kinda, qs_kindb
3260 : TYPE(xtb_atom_type), POINTER :: xtb_parameter_a, xtb_parameter_b
3261 :
3262 1232 : NULLIFY (logger)
3263 1232 : logger => cp_get_default_logger()
3264 1232 : IF (BTEST(cp_print_key_should_output(logger%iter_info, subsys_section, &
3265 : "PRINT%KINDS/POTENTIAL"), cp_p_file)) THEN
3266 :
3267 0 : io_unit = cp_print_key_unit_nr(logger, subsys_section, "PRINT%KINDS", extension=".Log")
3268 0 : IF (io_unit > 0) THEN
3269 :
3270 0 : WRITE (io_unit, "(/,T2,A)") "xTB| Kab parameters"
3271 0 : nkind = SIZE(qs_kind_set)
3272 0 : DO ikind = 1, nkind
3273 0 : qs_kinda => qs_kind_set(ikind)
3274 0 : CALL get_qs_kind(qs_kinda, xtb_parameter=xtb_parameter_a)
3275 0 : CALL get_xtb_atom_param(xtb_parameter_a, aname=aname, z=za)
3276 0 : DO jkind = ikind, nkind
3277 0 : qs_kindb => qs_kind_set(jkind)
3278 0 : CALL get_qs_kind(qs_kindb, xtb_parameter=xtb_parameter_b)
3279 0 : CALL get_xtb_atom_param(xtb_parameter_b, aname=bname, z=zb)
3280 : WRITE (io_unit, "(A,T10,A15,T25,A15,T71,F10.3)") &
3281 0 : " Kab:", TRIM(aname), TRIM(bname), xtb_set_kab(za, zb, xtb_control)
3282 : END DO
3283 : END DO
3284 0 : WRITE (io_unit, *)
3285 :
3286 : END IF
3287 :
3288 0 : CALL cp_print_key_finished_output(io_unit, logger, subsys_section, "PRINT%KINDS")
3289 : END IF
3290 :
3291 1232 : END SUBROUTINE write_xtb_kab_param
3292 :
3293 : ! **************************************************************************************************
3294 : !> \brief Set the components of an atomic kind data set.
3295 : !> \param qs_kind ...
3296 : !> \param paw_atom ...
3297 : !> \param ghost ...
3298 : !> \param floating ...
3299 : !> \param hard_radius ...
3300 : !> \param hard0_radius ...
3301 : !> \param covalent_radius ...
3302 : !> \param vdw_radius ...
3303 : !> \param lmax_rho0 ...
3304 : !> \param zeff ...
3305 : !> \param no_optimize ...
3306 : !> \param dispersion ...
3307 : !> \param u_minus_j ...
3308 : !> \param hund_j ...
3309 : !> \param reltmat ...
3310 : !> \param dftb_parameter ...
3311 : !> \param xtb_parameter ...
3312 : !> \param elec_conf ...
3313 : !> \param pao_basis_size ...
3314 : ! **************************************************************************************************
3315 31597 : SUBROUTINE set_qs_kind(qs_kind, paw_atom, ghost, floating, hard_radius, hard0_radius, &
3316 : covalent_radius, vdw_radius, lmax_rho0, zeff, &
3317 : no_optimize, dispersion, u_minus_j, hund_j, reltmat, &
3318 : dftb_parameter, xtb_parameter, &
3319 31597 : elec_conf, pao_basis_size)
3320 :
3321 : TYPE(qs_kind_type), INTENT(INOUT) :: qs_kind
3322 : LOGICAL, INTENT(IN), OPTIONAL :: paw_atom, ghost, floating
3323 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: hard_radius, hard0_radius, &
3324 : covalent_radius, vdw_radius
3325 : INTEGER, INTENT(IN), OPTIONAL :: lmax_rho0
3326 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: zeff
3327 : LOGICAL, INTENT(IN), OPTIONAL :: no_optimize
3328 : TYPE(qs_atom_dispersion_type), OPTIONAL, POINTER :: dispersion
3329 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: u_minus_j, hund_j
3330 : REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: reltmat
3331 : TYPE(qs_dftb_atom_type), OPTIONAL, POINTER :: dftb_parameter
3332 : TYPE(xtb_atom_type), OPTIONAL, POINTER :: xtb_parameter
3333 : INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: elec_conf
3334 : INTEGER, INTENT(IN), OPTIONAL :: pao_basis_size
3335 :
3336 31597 : IF (PRESENT(dftb_parameter)) qs_kind%dftb_parameter => dftb_parameter
3337 31597 : IF (PRESENT(xtb_parameter)) qs_kind%xtb_parameter => xtb_parameter
3338 31597 : IF (PRESENT(elec_conf)) THEN
3339 17159 : IF (ASSOCIATED(qs_kind%elec_conf)) THEN
3340 0 : DEALLOCATE (qs_kind%elec_conf)
3341 : END IF
3342 51477 : ALLOCATE (qs_kind%elec_conf(0:SIZE(elec_conf) - 1))
3343 62549 : qs_kind%elec_conf(:) = elec_conf(:)
3344 : END IF
3345 31597 : IF (PRESENT(paw_atom)) qs_kind%paw_atom = paw_atom
3346 31597 : IF (PRESENT(hard_radius)) qs_kind%hard_radius = hard_radius
3347 31597 : IF (PRESENT(hard0_radius)) qs_kind%hard0_radius = hard0_radius
3348 31597 : IF (PRESENT(covalent_radius)) qs_kind%covalent_radius = covalent_radius
3349 31597 : IF (PRESENT(vdw_radius)) qs_kind%vdw_radius = vdw_radius
3350 31597 : IF (PRESENT(lmax_rho0)) qs_kind%lmax_rho0 = lmax_rho0
3351 31597 : IF (PRESENT(zeff)) THEN
3352 0 : IF (ASSOCIATED(qs_kind%all_potential)) THEN
3353 0 : CALL set_potential(potential=qs_kind%all_potential, zeff=zeff)
3354 0 : ELSE IF (ASSOCIATED(qs_kind%gth_potential)) THEN
3355 0 : CALL set_potential(potential=qs_kind%gth_potential, zeff=zeff)
3356 0 : ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
3357 0 : CALL set_potential(potential=qs_kind%sgp_potential, zeff=zeff)
3358 0 : ELSE IF (ASSOCIATED(qs_kind%cneo_potential)) THEN
3359 0 : CPABORT("CNEO potential ZEFF should not be manually set")
3360 : END IF
3361 : END IF
3362 31597 : IF (PRESENT(ghost)) qs_kind%ghost = ghost
3363 :
3364 31597 : IF (PRESENT(floating)) qs_kind%floating = floating
3365 :
3366 31597 : IF (PRESENT(no_optimize)) qs_kind%no_optimize = no_optimize
3367 :
3368 31597 : IF (PRESENT(dispersion)) qs_kind%dispersion => dispersion
3369 :
3370 31597 : IF (PRESENT(u_minus_j)) THEN
3371 476 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
3372 476 : qs_kind%dft_plus_u%u_minus_j = u_minus_j
3373 : END IF
3374 : END IF
3375 :
3376 31597 : IF (PRESENT(hund_j)) THEN
3377 0 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
3378 0 : qs_kind%dft_plus_u%hund_j = hund_j
3379 : END IF
3380 : END IF
3381 :
3382 31597 : IF (PRESENT(reltmat)) qs_kind%reltmat => reltmat
3383 :
3384 31597 : IF (PRESENT(pao_basis_size)) qs_kind%pao_basis_size = pao_basis_size
3385 :
3386 31597 : END SUBROUTINE set_qs_kind
3387 :
3388 : ! **************************************************************************************************
3389 : !> \brief Write an atomic kind data set to the output unit.
3390 : !> \param qs_kind ...
3391 : !> \param kind_number ...
3392 : !> \param output_unit ...
3393 : !> \par History
3394 : !> Creation (09.02.2002,MK)
3395 : ! **************************************************************************************************
3396 4042 : SUBROUTINE write_qs_kind(qs_kind, kind_number, output_unit)
3397 :
3398 : TYPE(qs_kind_type), POINTER :: qs_kind
3399 : INTEGER, INTENT(in) :: kind_number, output_unit
3400 :
3401 : CHARACTER(LEN=3) :: yon
3402 : CHARACTER(LEN=default_string_length) :: basis_type, bstring
3403 : INTEGER :: ibas
3404 : LOGICAL :: do_print
3405 : TYPE(gto_basis_set_type), POINTER :: tmp_basis
3406 :
3407 4042 : IF (output_unit > 0) THEN
3408 :
3409 4042 : IF (ASSOCIATED(qs_kind)) THEN
3410 : WRITE (UNIT=output_unit, FMT="(/,T2,I2,A,T57,A,T75,I6)") &
3411 4042 : kind_number, ". Atomic kind: "//TRIM(qs_kind%name), &
3412 8084 : "Number of atoms: ", qs_kind%natom
3413 :
3414 84882 : DO ibas = 1, SIZE(qs_kind%basis_sets, 1)
3415 80840 : NULLIFY (tmp_basis)
3416 : CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis, &
3417 80840 : inumbas=ibas, basis_type=basis_type)
3418 80840 : do_print = .TRUE.
3419 75487 : SELECT CASE (basis_type)
3420 : CASE DEFAULT
3421 75487 : bstring = "Basis Set"
3422 3937 : do_print = .FALSE.
3423 : CASE ("ORB")
3424 3937 : bstring = "Orbital Basis Set"
3425 : CASE ("ORB_SOFT")
3426 629 : bstring = "GAPW Soft Basis Set"
3427 0 : do_print = .FALSE.
3428 : CASE ("AUX")
3429 0 : bstring = "Auxiliary Basis Set"
3430 : CASE ("MIN")
3431 0 : bstring = "Minimal Basis Set"
3432 : CASE ("RI_AUX")
3433 420 : bstring = "RI Auxiliary Basis Set"
3434 : CASE ("AUX_FIT")
3435 250 : bstring = "Auxiliary Fit Basis Set"
3436 : CASE ("LRI_AUX")
3437 19 : bstring = "LRI Basis Set"
3438 : CASE ("P_LRI_AUX")
3439 4 : bstring = "LRI Basis Set for TDDFPT"
3440 : CASE ("RI_XAS")
3441 0 : bstring = "RI XAS Basis Set"
3442 : CASE ("RI_HFX")
3443 86 : bstring = "RI HFX Basis Set"
3444 : CASE ("NUC")
3445 4 : bstring = "Nuclear Basis Set"
3446 4 : do_print = .FALSE.
3447 : CASE ("NUC_SOFT")
3448 4 : bstring = "Nuclear Soft Basis Set"
3449 80840 : do_print = .FALSE.
3450 : END SELECT
3451 :
3452 4042 : IF (do_print) THEN
3453 4716 : CALL write_orb_basis_set(tmp_basis, output_unit, bstring)
3454 : END IF
3455 :
3456 : END DO
3457 :
3458 4042 : IF (qs_kind%ghost) THEN
3459 : WRITE (UNIT=output_unit, FMT="(/,T6,A)") &
3460 11 : "The atoms of this atomic kind are GHOST atoms!"
3461 : END IF
3462 4042 : IF (qs_kind%monovalent) THEN
3463 : WRITE (UNIT=output_unit, FMT="(/,T6,A)") &
3464 1 : "The atoms of this atomic kind are MONOVALENT!"
3465 : END IF
3466 4042 : IF (qs_kind%floating) THEN
3467 : WRITE (UNIT=output_unit, FMT="(/,T6,A)") &
3468 0 : "The atoms of this atomic kind are FLOATING BASIS FUNCTIONS."
3469 : END IF
3470 4042 : IF (qs_kind%covalent_radius > 0.0_dp) THEN
3471 : WRITE (UNIT=output_unit, FMT="(/,T8,A,T71,F10.3)") &
3472 2734 : "Atomic covalent radius [Angstrom]:", &
3473 5468 : qs_kind%covalent_radius*angstrom
3474 : END IF
3475 4042 : IF (qs_kind%vdw_radius > 0.0_dp) THEN
3476 : WRITE (UNIT=output_unit, FMT="(/,T8,A,T71,F10.3)") &
3477 2734 : "Atomic van der Waals radius [Angstrom]:", &
3478 5468 : qs_kind%vdw_radius*angstrom
3479 : END IF
3480 4042 : IF (qs_kind%paw_atom) THEN
3481 : WRITE (UNIT=output_unit, FMT="(/,T6,A)") &
3482 507 : "The atoms of this atomic kind are PAW atoms (GAPW):"
3483 : WRITE (UNIT=output_unit, FMT="(T8,A,T71,F10.3)") &
3484 507 : "Hard Gaussian function radius:", qs_kind%hard_radius, &
3485 507 : "Rho0 radius:", qs_kind%hard0_radius, &
3486 507 : "Maximum GTO radius used for PAW projector construction:", &
3487 1014 : qs_kind%max_rad_local
3488 507 : NULLIFY (tmp_basis)
3489 : CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis, &
3490 507 : basis_type="ORB_SOFT")
3491 507 : CALL write_orb_basis_set(tmp_basis, output_unit, "GAPW Soft Basis Set")
3492 : END IF
3493 : ! Potentials
3494 4042 : IF (ASSOCIATED(qs_kind%all_potential)) CALL write_potential(qs_kind%all_potential, output_unit)
3495 4042 : IF (ASSOCIATED(qs_kind%gth_potential)) CALL write_potential(qs_kind%gth_potential, output_unit)
3496 4042 : IF (ASSOCIATED(qs_kind%sgp_potential)) CALL write_potential(qs_kind%sgp_potential, output_unit)
3497 4042 : IF (ASSOCIATED(qs_kind%tnadd_potential)) CALL write_potential(qs_kind%tnadd_potential, output_unit)
3498 4042 : IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
3499 : WRITE (UNIT=output_unit, FMT="(/,T6,A,/,T8,A,T76,I5,/,T8,A,T73,F8.3,/,T8,A,T73,F8.3)") &
3500 19 : "A DFT+U correction is applied to atoms of this atomic kind:", &
3501 19 : "Angular quantum momentum number L:", qs_kind%dft_plus_u%l, &
3502 19 : "U(eff) = (U - J) value in [eV]:", qs_kind%dft_plus_u%u_minus_j_target*evolt, &
3503 38 : "Hund J value in [eV]:", qs_kind%dft_plus_u%hund_j*evolt
3504 19 : IF (qs_kind%dft_plus_u%u_ramping > 0.0_dp) THEN
3505 4 : IF (qs_kind%dft_plus_u%init_u_ramping_each_scf) THEN
3506 2 : yon = "YES"
3507 : ELSE
3508 2 : yon = " NO"
3509 : END IF
3510 : WRITE (UNIT=output_unit, FMT="(T8,A,T73,F8.3,/,T8,A,T73,ES8.1,/,T8,A,T78,A3)") &
3511 4 : "Increment for U ramping in [eV]:", qs_kind%dft_plus_u%u_ramping*evolt, &
3512 4 : "SCF threshold value for U ramping:", qs_kind%dft_plus_u%eps_u_ramping, &
3513 8 : "Set U ramping value to zero before each wavefunction optimisation:", yon
3514 : END IF
3515 19 : IF (ASSOCIATED(qs_kind%dft_plus_u%orbitals)) THEN
3516 : WRITE (UNIT=output_unit, FMT="(T8,A)") &
3517 2 : "An initial orbital occupation is requested:"
3518 2 : IF (ASSOCIATED(qs_kind%dft_plus_u%nelec)) THEN
3519 4 : IF (ANY(qs_kind%dft_plus_u%nelec(:) >= 0.5_dp)) THEN
3520 0 : IF (SIZE(qs_kind%dft_plus_u%nelec) > 1) THEN
3521 : WRITE (UNIT=output_unit, FMT="(T9,A,T75,F6.2)") &
3522 0 : "Number of alpha electrons:", &
3523 0 : qs_kind%dft_plus_u%nelec(1), &
3524 0 : "Number of beta electrons:", &
3525 0 : qs_kind%dft_plus_u%nelec(2)
3526 : ELSE
3527 : WRITE (UNIT=output_unit, FMT="(T9,A,T75,F6.2)") &
3528 0 : "Number of electrons:", &
3529 0 : qs_kind%dft_plus_u%nelec(1)
3530 : END IF
3531 : END IF
3532 : END IF
3533 : WRITE (UNIT=output_unit, FMT="(T9,A,(T78,I3))") &
3534 2 : "Preferred (initial) orbital occupation order (orbital M values):", &
3535 10 : qs_kind%dft_plus_u%orbitals(:)
3536 : WRITE (UNIT=output_unit, FMT="(T9,A,T71,ES10.3,/,T9,A,T76,I5)") &
3537 2 : "Threshold value for the SCF convergence criterion:", &
3538 2 : qs_kind%dft_plus_u%eps_scf, &
3539 2 : "Number of initial SCF iterations:", &
3540 4 : qs_kind%dft_plus_u%max_scf
3541 2 : IF (qs_kind%dft_plus_u%smear) THEN
3542 : WRITE (UNIT=output_unit, FMT="(T9,A)") &
3543 2 : "A smearing of the orbital occupations will be performed"
3544 : END IF
3545 : END IF
3546 : END IF
3547 4042 : IF (ASSOCIATED(qs_kind%cneo_potential)) THEN
3548 : WRITE (UNIT=output_unit, FMT="(/,T6,A)") &
3549 4 : "The nuclei of this atomic kind are quantum mechanical (CNEO)"
3550 4 : CALL write_cneo_potential(qs_kind%cneo_potential, output_unit)
3551 4 : NULLIFY (tmp_basis)
3552 : CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis, &
3553 4 : basis_type="NUC")
3554 4 : CALL write_orb_basis_set(tmp_basis, output_unit, "Nuclear Basis Set")
3555 4 : NULLIFY (tmp_basis)
3556 : CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis, &
3557 4 : basis_type="NUC_SOFT")
3558 4 : CALL write_orb_basis_set(tmp_basis, output_unit, "Nuclear Soft Basis Set")
3559 : END IF
3560 : ELSE
3561 0 : CPABORT("write_qs_kind requires qs_kind to be associated")
3562 : END IF
3563 :
3564 : END IF
3565 :
3566 4042 : END SUBROUTINE write_qs_kind
3567 :
3568 : ! **************************************************************************************************
3569 : !> \brief Write an atomic kind set data set to the output unit.
3570 : !> \param qs_kind_set ...
3571 : !> \param subsys_section ...
3572 : !> \par History
3573 : !> Creation (09.02.2002,MK)
3574 : ! **************************************************************************************************
3575 9320 : SUBROUTINE write_qs_kind_set(qs_kind_set, subsys_section)
3576 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3577 : TYPE(section_vals_type), POINTER :: subsys_section
3578 :
3579 : CHARACTER(len=*), PARAMETER :: routineN = 'write_qs_kind_set'
3580 :
3581 : INTEGER :: handle, ikind, nkind, output_unit
3582 : TYPE(cp_logger_type), POINTER :: logger
3583 : TYPE(qs_kind_type), POINTER :: qs_kind
3584 :
3585 9320 : CALL timeset(routineN, handle)
3586 :
3587 9320 : NULLIFY (logger)
3588 9320 : logger => cp_get_default_logger()
3589 : output_unit = cp_print_key_unit_nr(logger, subsys_section, &
3590 9320 : "PRINT%KINDS", extension=".Log")
3591 9320 : IF (output_unit > 0) THEN
3592 2194 : IF (ASSOCIATED(qs_kind_set)) THEN
3593 2194 : WRITE (UNIT=output_unit, FMT="(/,/,T2,A)") "ATOMIC KIND INFORMATION"
3594 2194 : nkind = SIZE(qs_kind_set)
3595 6236 : DO ikind = 1, nkind
3596 4042 : qs_kind => qs_kind_set(ikind)
3597 6236 : CALL write_qs_kind(qs_kind, ikind, output_unit)
3598 : END DO
3599 : ELSE
3600 0 : CPABORT("write_qs_kind_set requires qs_kind_set to be associated")
3601 : END IF
3602 : END IF
3603 :
3604 : CALL cp_print_key_finished_output(output_unit, logger, subsys_section, &
3605 9320 : "PRINT%KINDS")
3606 :
3607 9320 : CALL timestop(handle)
3608 :
3609 9320 : END SUBROUTINE write_qs_kind_set
3610 :
3611 : ! **************************************************************************************************
3612 : !> \brief Write all the GTO basis sets of an atomic kind set to the output
3613 : !> unit (for the printing of the unnormalized basis sets as read from
3614 : !> database).
3615 : !> \param qs_kind_set ...
3616 : !> \param subsys_section ...
3617 : !> \par History
3618 : !> Creation (17.01.2002,MK)
3619 : ! **************************************************************************************************
3620 9296 : SUBROUTINE write_gto_basis_sets(qs_kind_set, subsys_section)
3621 :
3622 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3623 : TYPE(section_vals_type), POINTER :: subsys_section
3624 :
3625 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_gto_basis_sets'
3626 :
3627 : CHARACTER(LEN=default_string_length) :: basis_type, bstring
3628 : INTEGER :: handle, ibas, ikind, nkind, output_unit
3629 : TYPE(cp_logger_type), POINTER :: logger
3630 : TYPE(gto_basis_set_type), POINTER :: tmp_basis
3631 : TYPE(qs_kind_type), POINTER :: qs_kind
3632 :
3633 9296 : CALL timeset(routineN, handle)
3634 :
3635 9296 : NULLIFY (logger)
3636 9296 : logger => cp_get_default_logger()
3637 : output_unit = cp_print_key_unit_nr(logger, subsys_section, &
3638 : "PRINT%KINDS/BASIS_SET", &
3639 9296 : extension=".Log")
3640 9296 : IF (output_unit > 0) THEN
3641 63 : IF (ASSOCIATED(qs_kind_set)) THEN
3642 : WRITE (UNIT=output_unit, FMT="(/,/,T2,A)") &
3643 63 : "BASIS SET INFORMATION (Unnormalised Gaussian-type functions)"
3644 63 : nkind = SIZE(qs_kind_set)
3645 181 : DO ikind = 1, nkind
3646 118 : qs_kind => qs_kind_set(ikind)
3647 : WRITE (UNIT=output_unit, FMT="(/,T2,I2,A)") &
3648 118 : ikind, ". Atomic kind: "//TRIM(qs_kind%name)
3649 :
3650 2541 : DO ibas = 1, SIZE(qs_kind%basis_sets, 1)
3651 2360 : NULLIFY (tmp_basis)
3652 : CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis, &
3653 2360 : inumbas=ibas, basis_type=basis_type)
3654 2360 : IF (basis_type == "") CYCLE
3655 11 : SELECT CASE (basis_type)
3656 : CASE DEFAULT
3657 11 : bstring = "Basis Set"
3658 : CASE ("ORB")
3659 118 : bstring = "Orbital Basis Set"
3660 : CASE ("ORB_SOFT")
3661 11 : bstring = "GAPW Soft Basis Set"
3662 : CASE ("AUX")
3663 0 : bstring = "Auxiliary Basis Set"
3664 : CASE ("MIN")
3665 0 : bstring = "Minimal Basis Set"
3666 : CASE ("RI_AUX")
3667 0 : bstring = "RI Auxiliary Basis Set"
3668 : CASE ("AUX_FIT")
3669 0 : bstring = "Auxiliary Fit Basis Set"
3670 : CASE ("LRI_AUX")
3671 2 : bstring = "LRI Basis Set"
3672 : CASE ("P_LRI_AUX")
3673 0 : bstring = "LRI Basis Set for TDDFPT"
3674 : CASE ("RI_HFX")
3675 0 : bstring = "RI HFX Basis Set"
3676 : CASE ("NUC")
3677 0 : bstring = "Nuclear Basis Set"
3678 0 : IF (.NOT. ASSOCIATED(qs_kind%cneo_potential)) NULLIFY (tmp_basis)
3679 : CASE ("NUC_SOFT")
3680 0 : bstring = "Nuclear Soft Basis Set"
3681 142 : IF (.NOT. ASSOCIATED(qs_kind%cneo_potential)) NULLIFY (tmp_basis)
3682 : END SELECT
3683 :
3684 260 : IF (ASSOCIATED(tmp_basis)) CALL write_gto_basis_set(tmp_basis, output_unit, bstring)
3685 :
3686 : END DO
3687 :
3688 : END DO
3689 : ELSE
3690 0 : CPABORT("write_gto_basis_sets require qs_kind_set to be associated")
3691 : END IF
3692 : END IF
3693 :
3694 : CALL cp_print_key_finished_output(output_unit, logger, subsys_section, &
3695 9296 : "PRINT%KINDS/BASIS_SET")
3696 :
3697 9296 : CALL timestop(handle)
3698 :
3699 9296 : END SUBROUTINE write_gto_basis_sets
3700 :
3701 : ! **************************************************************************************************
3702 : !> \brief ...
3703 : !> \param atomic_kind ...
3704 : !> \param qs_kind ...
3705 : !> \param ncalc ...
3706 : !> \param ncore ...
3707 : !> \param nelem ...
3708 : !> \param edelta ...
3709 : !> \param rks ...
3710 : ! **************************************************************************************************
3711 110494 : SUBROUTINE init_atom_electronic_state(atomic_kind, qs_kind, ncalc, ncore, nelem, edelta, rks)
3712 :
3713 : TYPE(atomic_kind_type), INTENT(IN) :: atomic_kind
3714 : TYPE(qs_kind_type), INTENT(IN) :: qs_kind
3715 : INTEGER, DIMENSION(0:lmat, 10), INTENT(OUT) :: ncalc, ncore, nelem
3716 : REAL(KIND=dp), DIMENSION(0:lmat, 10, 2), &
3717 : INTENT(OUT) :: edelta
3718 : LOGICAL, INTENT(IN), OPTIONAL :: rks
3719 :
3720 : INTEGER :: i, ii, is, l, ll, ne, nn, z
3721 55247 : INTEGER, DIMENSION(:), POINTER :: econf
3722 55247 : INTEGER, DIMENSION(:, :), POINTER :: addel, laddel, naddel
3723 : LOGICAL :: bs_occupation, monovalent
3724 : REAL(KIND=dp) :: dmag, magnetization
3725 : TYPE(gth_potential_type), POINTER :: gth_potential
3726 : TYPE(sgp_potential_type), POINTER :: sgp_potential
3727 :
3728 55247 : CALL get_atomic_kind(atomic_kind, z=z)
3729 55247 : NULLIFY (gth_potential)
3730 : CALL get_qs_kind(qs_kind, &
3731 : gth_potential=gth_potential, &
3732 : sgp_potential=sgp_potential, &
3733 : magnetization=magnetization, &
3734 : bs_occupation=bs_occupation, &
3735 : monovalent=monovalent, &
3736 55247 : addel=addel, laddel=laddel, naddel=naddel)
3737 :
3738 55247 : IF (PRESENT(rks)) THEN
3739 10970 : IF (rks) THEN
3740 56 : magnetization = 0.0_dp
3741 56 : bs_occupation = .FALSE.
3742 : END IF
3743 : END IF
3744 :
3745 : ! electronic state
3746 55247 : nelem = 0
3747 55247 : ncore = 0
3748 55247 : ncalc = 0
3749 55247 : edelta = 0.0_dp
3750 55247 : IF (monovalent) THEN
3751 4 : ncalc(0, 1) = 1
3752 4 : nelem(0, 1) = 1
3753 55243 : ELSE IF (ASSOCIATED(gth_potential)) THEN
3754 31037 : CALL get_potential(gth_potential, elec_conf=econf)
3755 31037 : CALL set_pseudo_state(econf, z, ncalc, ncore, nelem)
3756 24206 : ELSE IF (ASSOCIATED(sgp_potential)) THEN
3757 228 : CALL get_potential(sgp_potential, elec_conf=econf)
3758 228 : CALL set_pseudo_state(econf, z, ncalc, ncore, nelem)
3759 : ELSE
3760 119890 : DO l = 0, MIN(lmat, UBOUND(ptable(z)%e_conv, 1))
3761 95912 : ll = 2*(2*l + 1)
3762 95912 : nn = ptable(z)%e_conv(l)
3763 95912 : ii = 0
3764 23978 : DO
3765 135230 : ii = ii + 1
3766 135230 : IF (nn <= ll) THEN
3767 95912 : nelem(l, ii) = nn
3768 : EXIT
3769 : ELSE
3770 39318 : nelem(l, ii) = ll
3771 39318 : nn = nn - ll
3772 : END IF
3773 : END DO
3774 : END DO
3775 1702438 : ncalc = nelem - ncore
3776 : END IF
3777 :
3778 : ! readjust the occupation number of the orbitals as requested by user
3779 : ! this is done to break symmetry (bs) and bias the initial guess
3780 : ! to the pre-defined multiplicity/charge state of the atom
3781 55247 : IF (bs_occupation) THEN
3782 672 : DO is = 1, 2
3783 1216 : DO i = 1, SIZE(addel, 1)
3784 544 : ne = addel(i, is)
3785 544 : l = laddel(i, is)
3786 544 : nn = naddel(i, is) - l
3787 992 : IF (ne /= 0) THEN
3788 516 : IF (nn == 0) THEN
3789 0 : DO ii = SIZE(nelem, 2), 1, -1
3790 0 : IF (ncalc(l, ii) > 0) THEN
3791 0 : IF ((ncalc(l, ii) + ne) < 2*(2*l + 1) + 1) THEN
3792 0 : edelta(l, ii, is) = edelta(l, ii, is) + ne
3793 0 : nn = ii
3794 : ELSE
3795 0 : edelta(l, ii + 1, is) = edelta(l, ii + 1, is) + ne
3796 0 : nn = ii + 1
3797 : END IF
3798 : EXIT
3799 0 : ELSE IF (ii == 1) THEN
3800 0 : edelta(l, ii, is) = edelta(l, ii, is) + ne
3801 0 : nn = ii
3802 : END IF
3803 : END DO
3804 : ELSE
3805 516 : edelta(l, nn, is) = edelta(l, nn, is) + ne
3806 : END IF
3807 516 : IF (ncalc(l, nn) + edelta(l, nn, is) < 0) THEN
3808 0 : edelta(l, nn, is) = -ncalc(l, nn)
3809 : END IF
3810 : END IF
3811 : END DO
3812 : END DO
3813 32032 : edelta = 0.5_dp*edelta
3814 55023 : ELSE IF (magnetization /= 0.0_dp) THEN
3815 40 : dmag = 0.5_dp*ABS(magnetization)
3816 94 : DO l = 0, MIN(lmat, UBOUND(ptable(z)%e_conv, 1))
3817 94 : ll = 2*(2*l + 1)
3818 94 : ii = 0
3819 660 : DO i = 1, SIZE(ncalc, 2)
3820 606 : IF (ncalc(l, i) == 0) CYCLE
3821 116 : IF (ncalc(l, i) == ll) CYCLE
3822 100 : IF (ncalc(l, i) > dmag .AND. (ll - ncalc(l, i)) > dmag) THEN
3823 : ii = i
3824 : EXIT
3825 : END IF
3826 : END DO
3827 94 : IF (ii /= 0) THEN
3828 40 : edelta(l, ii, 1) = magnetization*0.5_dp
3829 40 : edelta(l, ii, 2) = -magnetization*0.5_dp
3830 40 : EXIT
3831 : END IF
3832 : END DO
3833 40 : IF (ii == 0) THEN
3834 : CALL cp_abort(__LOCATION__, &
3835 0 : "Magnetization value cannot be imposed for this atom type")
3836 : END IF
3837 : END IF
3838 :
3839 55247 : IF (qs_kind%ghost .OR. qs_kind%floating) THEN
3840 408 : nelem = 0
3841 408 : ncore = 0
3842 408 : ncalc = 0
3843 408 : edelta = 0.0_dp
3844 : END IF
3845 :
3846 55247 : END SUBROUTINE init_atom_electronic_state
3847 :
3848 : ! **************************************************************************************************
3849 : !> \brief ...
3850 : !> \param econf ...
3851 : !> \param z ...
3852 : !> \param ncalc ...
3853 : !> \param ncore ...
3854 : !> \param nelem ...
3855 : ! **************************************************************************************************
3856 31325 : SUBROUTINE set_pseudo_state(econf, z, ncalc, ncore, nelem)
3857 : INTEGER, DIMENSION(:), POINTER :: econf
3858 : INTEGER, INTENT(IN) :: z
3859 : INTEGER, DIMENSION(0:lmat, 10), INTENT(OUT) :: ncalc, ncore, nelem
3860 :
3861 : CHARACTER(LEN=default_string_length) :: message
3862 : INTEGER :: ii, iounit, l, ll, lmin, nc, nn
3863 : INTEGER, DIMENSION(0:lmat) :: econfx
3864 : TYPE(cp_logger_type), POINTER :: logger
3865 :
3866 31325 : NULLIFY (logger)
3867 31325 : logger => cp_get_default_logger()
3868 31325 : iounit = cp_logger_get_default_io_unit(logger)
3869 :
3870 31325 : econfx = 0
3871 89026 : econfx(0:SIZE(econf) - 1) = econf
3872 89026 : IF (SUM(econf) >= 0) THEN
3873 88958 : lmin = MIN(lmat, UBOUND(ptable(z)%e_conv, 1))
3874 : ! number of core electrons
3875 88958 : nc = z - SUM(econf)
3876 : ! setup ncore
3877 31291 : ncore = 0
3878 11284 : SELECT CASE (nc)
3879 : CASE (0)
3880 : CASE (2)
3881 11284 : ncore(0, 1) = 2
3882 : CASE (10)
3883 2744 : ncore(0, 1) = 2
3884 2744 : ncore(0, 2) = 2
3885 2744 : ncore(1, 1) = 6
3886 : CASE (18)
3887 58 : ncore(0, 1) = 2
3888 58 : ncore(0, 2) = 2
3889 58 : ncore(0, 3) = 2
3890 58 : ncore(1, 1) = 6
3891 58 : ncore(1, 2) = 6
3892 : CASE (28)
3893 66 : ncore(0, 1) = 2
3894 66 : ncore(0, 2) = 2
3895 66 : ncore(0, 3) = 2
3896 66 : ncore(1, 1) = 6
3897 66 : ncore(1, 2) = 6
3898 66 : ncore(2, 1) = 10
3899 : CASE (36)
3900 0 : ncore(0, 1) = 2
3901 0 : ncore(0, 2) = 2
3902 0 : ncore(0, 3) = 2
3903 0 : ncore(0, 4) = 2
3904 0 : ncore(1, 1) = 6
3905 0 : ncore(1, 2) = 6
3906 0 : ncore(1, 3) = 6
3907 0 : ncore(2, 1) = 10
3908 : CASE (46)
3909 92 : ncore(0, 1) = 2
3910 92 : ncore(0, 2) = 2
3911 92 : ncore(0, 3) = 2
3912 92 : ncore(0, 4) = 2
3913 92 : ncore(1, 1) = 6
3914 92 : ncore(1, 2) = 6
3915 92 : ncore(1, 3) = 6
3916 92 : ncore(2, 1) = 10
3917 92 : ncore(2, 2) = 10
3918 : CASE (54)
3919 4 : ncore(0, 1) = 2
3920 4 : ncore(0, 2) = 2
3921 4 : ncore(0, 3) = 2
3922 4 : ncore(0, 4) = 2
3923 4 : ncore(0, 5) = 2
3924 4 : ncore(1, 1) = 6
3925 4 : ncore(1, 2) = 6
3926 4 : ncore(1, 3) = 6
3927 4 : ncore(1, 4) = 6
3928 4 : ncore(2, 1) = 10
3929 4 : ncore(2, 2) = 10
3930 : CASE (60)
3931 34 : ncore(0, 1) = 2
3932 34 : ncore(0, 2) = 2
3933 34 : ncore(0, 3) = 2
3934 34 : ncore(0, 4) = 2
3935 34 : ncore(1, 1) = 6
3936 34 : ncore(1, 2) = 6
3937 34 : ncore(1, 3) = 6
3938 34 : ncore(2, 1) = 10
3939 34 : ncore(2, 2) = 10
3940 34 : ncore(3, 1) = 14
3941 : CASE (68)
3942 250 : ncore(0, 1) = 2
3943 250 : ncore(0, 2) = 2
3944 250 : ncore(0, 3) = 2
3945 250 : ncore(0, 4) = 2
3946 250 : ncore(0, 5) = 2
3947 250 : ncore(1, 1) = 6
3948 250 : ncore(1, 2) = 6
3949 250 : ncore(1, 3) = 6
3950 250 : ncore(1, 4) = 6
3951 250 : ncore(2, 1) = 10
3952 250 : ncore(2, 2) = 10
3953 250 : ncore(3, 1) = 14
3954 : CASE (78)
3955 28 : ncore(0, 1) = 2
3956 28 : ncore(0, 2) = 2
3957 28 : ncore(0, 3) = 2
3958 28 : ncore(0, 4) = 2
3959 28 : ncore(0, 5) = 2
3960 28 : ncore(1, 1) = 6
3961 28 : ncore(1, 2) = 6
3962 28 : ncore(1, 3) = 6
3963 28 : ncore(1, 4) = 6
3964 28 : ncore(2, 1) = 10
3965 28 : ncore(2, 2) = 10
3966 28 : ncore(2, 3) = 10
3967 28 : ncore(3, 1) = 14
3968 : ! 79 - 92 5f incore PP
3969 : CASE (79)
3970 0 : ncore(0, 1) = 2
3971 0 : ncore(0, 2) = 2
3972 0 : ncore(0, 3) = 2
3973 0 : ncore(0, 4) = 2
3974 0 : ncore(0, 5) = 2
3975 0 : ncore(1, 1) = 6
3976 0 : ncore(1, 2) = 6
3977 0 : ncore(1, 3) = 6
3978 0 : ncore(1, 4) = 6
3979 0 : ncore(2, 1) = 10
3980 0 : ncore(2, 2) = 10
3981 0 : ncore(2, 3) = 10
3982 0 : ncore(3, 1) = 14
3983 0 : ncore(3, 2) = 1
3984 : CASE (80)
3985 0 : ncore(0, 1) = 2
3986 0 : ncore(0, 2) = 2
3987 0 : ncore(0, 3) = 2
3988 0 : ncore(0, 4) = 2
3989 0 : ncore(0, 5) = 2
3990 0 : ncore(1, 1) = 6
3991 0 : ncore(1, 2) = 6
3992 0 : ncore(1, 3) = 6
3993 0 : ncore(1, 4) = 6
3994 0 : ncore(2, 1) = 10
3995 0 : ncore(2, 2) = 10
3996 0 : ncore(2, 3) = 10
3997 0 : ncore(3, 1) = 14
3998 0 : ncore(3, 2) = 2
3999 : CASE (81)
4000 0 : ncore(0, 1) = 2
4001 0 : ncore(0, 2) = 2
4002 0 : ncore(0, 3) = 2
4003 0 : ncore(0, 4) = 2
4004 0 : ncore(0, 5) = 2
4005 0 : ncore(1, 1) = 6
4006 0 : ncore(1, 2) = 6
4007 0 : ncore(1, 3) = 6
4008 0 : ncore(1, 4) = 6
4009 0 : ncore(2, 1) = 10
4010 0 : ncore(2, 2) = 10
4011 0 : ncore(2, 3) = 10
4012 0 : ncore(3, 1) = 14
4013 0 : ncore(3, 2) = 3
4014 : CASE (82)
4015 0 : ncore(0, 1) = 2
4016 0 : ncore(0, 2) = 2
4017 0 : ncore(0, 3) = 2
4018 0 : ncore(0, 4) = 2
4019 0 : ncore(0, 5) = 2
4020 0 : ncore(1, 1) = 6
4021 0 : ncore(1, 2) = 6
4022 0 : ncore(1, 3) = 6
4023 0 : ncore(1, 4) = 6
4024 0 : ncore(2, 1) = 10
4025 0 : ncore(2, 2) = 10
4026 0 : ncore(2, 3) = 10
4027 0 : ncore(3, 1) = 14
4028 0 : ncore(3, 2) = 4
4029 : CASE (83)
4030 0 : ncore(0, 1) = 2
4031 0 : ncore(0, 2) = 2
4032 0 : ncore(0, 3) = 2
4033 0 : ncore(0, 4) = 2
4034 0 : ncore(0, 5) = 2
4035 0 : ncore(1, 1) = 6
4036 0 : ncore(1, 2) = 6
4037 0 : ncore(1, 3) = 6
4038 0 : ncore(1, 4) = 6
4039 0 : ncore(2, 1) = 10
4040 0 : ncore(2, 2) = 10
4041 0 : ncore(2, 3) = 10
4042 0 : ncore(3, 1) = 14
4043 0 : ncore(3, 2) = 5
4044 : CASE (84)
4045 0 : ncore(0, 1) = 2
4046 0 : ncore(0, 2) = 2
4047 0 : ncore(0, 3) = 2
4048 0 : ncore(0, 4) = 2
4049 0 : ncore(0, 5) = 2
4050 0 : ncore(1, 1) = 6
4051 0 : ncore(1, 2) = 6
4052 0 : ncore(1, 3) = 6
4053 0 : ncore(1, 4) = 6
4054 0 : ncore(2, 1) = 10
4055 0 : ncore(2, 2) = 10
4056 0 : ncore(2, 3) = 10
4057 0 : ncore(3, 1) = 14
4058 0 : ncore(3, 2) = 6
4059 : CASE (85)
4060 0 : ncore(0, 1) = 2
4061 0 : ncore(0, 2) = 2
4062 0 : ncore(0, 3) = 2
4063 0 : ncore(0, 4) = 2
4064 0 : ncore(0, 5) = 2
4065 0 : ncore(1, 1) = 6
4066 0 : ncore(1, 2) = 6
4067 0 : ncore(1, 3) = 6
4068 0 : ncore(1, 4) = 6
4069 0 : ncore(2, 1) = 10
4070 0 : ncore(2, 2) = 10
4071 0 : ncore(2, 3) = 10
4072 0 : ncore(3, 1) = 14
4073 0 : ncore(3, 2) = 7
4074 : CASE (86)
4075 : ! this is not Rn core, add double assignment below
4076 0 : ncore(0, 1) = 2
4077 0 : ncore(0, 2) = 2
4078 0 : ncore(0, 3) = 2
4079 0 : ncore(0, 4) = 2
4080 0 : ncore(0, 5) = 2
4081 0 : ncore(1, 1) = 6
4082 0 : ncore(1, 2) = 6
4083 0 : ncore(1, 3) = 6
4084 0 : ncore(1, 4) = 6
4085 0 : ncore(2, 1) = 10
4086 0 : ncore(2, 2) = 10
4087 0 : ncore(2, 3) = 10
4088 0 : ncore(3, 1) = 14
4089 0 : ncore(3, 2) = 8
4090 : CASE (87)
4091 0 : ncore(0, 1) = 2
4092 0 : ncore(0, 2) = 2
4093 0 : ncore(0, 3) = 2
4094 0 : ncore(0, 4) = 2
4095 0 : ncore(0, 5) = 2
4096 0 : ncore(1, 1) = 6
4097 0 : ncore(1, 2) = 6
4098 0 : ncore(1, 3) = 6
4099 0 : ncore(1, 4) = 6
4100 0 : ncore(2, 1) = 10
4101 0 : ncore(2, 2) = 10
4102 0 : ncore(2, 3) = 10
4103 0 : ncore(3, 1) = 14
4104 0 : ncore(3, 2) = 9
4105 : CASE (88)
4106 0 : ncore(0, 1) = 2
4107 0 : ncore(0, 2) = 2
4108 0 : ncore(0, 3) = 2
4109 0 : ncore(0, 4) = 2
4110 0 : ncore(0, 5) = 2
4111 0 : ncore(1, 1) = 6
4112 0 : ncore(1, 2) = 6
4113 0 : ncore(1, 3) = 6
4114 0 : ncore(1, 4) = 6
4115 0 : ncore(2, 1) = 10
4116 0 : ncore(2, 2) = 10
4117 0 : ncore(2, 3) = 10
4118 0 : ncore(3, 1) = 14
4119 0 : ncore(3, 2) = 10
4120 : CASE (89)
4121 0 : ncore(0, 1) = 2
4122 0 : ncore(0, 2) = 2
4123 0 : ncore(0, 3) = 2
4124 0 : ncore(0, 4) = 2
4125 0 : ncore(0, 5) = 2
4126 0 : ncore(1, 1) = 6
4127 0 : ncore(1, 2) = 6
4128 0 : ncore(1, 3) = 6
4129 0 : ncore(1, 4) = 6
4130 0 : ncore(2, 1) = 10
4131 0 : ncore(2, 2) = 10
4132 0 : ncore(2, 3) = 10
4133 0 : ncore(3, 1) = 14
4134 0 : ncore(3, 2) = 11
4135 : CASE (90)
4136 0 : ncore(0, 1) = 2
4137 0 : ncore(0, 2) = 2
4138 0 : ncore(0, 3) = 2
4139 0 : ncore(0, 4) = 2
4140 0 : ncore(0, 5) = 2
4141 0 : ncore(1, 1) = 6
4142 0 : ncore(1, 2) = 6
4143 0 : ncore(1, 3) = 6
4144 0 : ncore(1, 4) = 6
4145 0 : ncore(2, 1) = 10
4146 0 : ncore(2, 2) = 10
4147 0 : ncore(2, 3) = 10
4148 0 : ncore(3, 1) = 14
4149 0 : ncore(3, 2) = 12
4150 : CASE (91)
4151 0 : ncore(0, 1) = 2
4152 0 : ncore(0, 2) = 2
4153 0 : ncore(0, 3) = 2
4154 0 : ncore(0, 4) = 2
4155 0 : ncore(0, 5) = 2
4156 0 : ncore(1, 1) = 6
4157 0 : ncore(1, 2) = 6
4158 0 : ncore(1, 3) = 6
4159 0 : ncore(1, 4) = 6
4160 0 : ncore(2, 1) = 10
4161 0 : ncore(2, 2) = 10
4162 0 : ncore(2, 3) = 10
4163 0 : ncore(3, 1) = 14
4164 0 : ncore(3, 2) = 13
4165 : CASE (92)
4166 0 : ncore(0, 1) = 2
4167 0 : ncore(0, 2) = 2
4168 0 : ncore(0, 3) = 2
4169 0 : ncore(0, 4) = 2
4170 0 : ncore(0, 5) = 2
4171 0 : ncore(1, 1) = 6
4172 0 : ncore(1, 2) = 6
4173 0 : ncore(1, 3) = 6
4174 0 : ncore(1, 4) = 6
4175 0 : ncore(2, 1) = 10
4176 0 : ncore(2, 2) = 10
4177 0 : ncore(2, 3) = 10
4178 0 : ncore(3, 1) = 14
4179 0 : ncore(3, 2) = 14
4180 : CASE DEFAULT
4181 31291 : ncore(0, 1) = -1
4182 : END SELECT
4183 : ! special cases of double assignments
4184 31291 : IF (z == 65 .AND. econfx(3) == 0) THEN
4185 : ! 4f in core for Tb
4186 4 : ncore = 0
4187 4 : ncore(0, 1) = -1
4188 : END IF
4189 : ! if there is still no core, check for special cases
4190 31291 : IF (ncore(0, 1) <= 0) THEN
4191 16735 : IF (z >= 58 .AND. z <= 71) THEN
4192 : ! 4f-in-core PPs for lanthanides
4193 280 : nc = z - SUM(econf)
4194 : ! setup ncore
4195 56 : ncore = 0
4196 0 : SELECT CASE (nc)
4197 : CASE (29:42)
4198 0 : ncore(0, 1) = 2
4199 0 : ncore(0, 2) = 2
4200 0 : ncore(0, 3) = 2
4201 0 : ncore(1, 1) = 6
4202 0 : ncore(1, 2) = 6
4203 0 : ncore(2, 1) = 10
4204 0 : ncore(3, 1) = nc - 28
4205 : message = "A small-core pseudopotential with 4f-in-core is used for the lanthanide "// &
4206 0 : TRIM(ptable(z)%symbol)
4207 0 : CPHINT(TRIM(message))
4208 : CASE (47:60)
4209 56 : ncore(0, 1) = 2
4210 56 : ncore(0, 2) = 2
4211 56 : ncore(0, 3) = 2
4212 56 : ncore(0, 4) = 2
4213 56 : ncore(1, 1) = 6
4214 56 : ncore(1, 2) = 6
4215 56 : ncore(1, 3) = 6
4216 56 : ncore(2, 1) = 10
4217 56 : ncore(2, 2) = 10
4218 56 : ncore(3, 1) = nc - 46
4219 : message = "A medium-core pseudopotential with 4f-in-core is used for the lanthanide "// &
4220 56 : TRIM(ptable(z)%symbol)
4221 56 : CPHINT(TRIM(message))
4222 : CASE DEFAULT
4223 56 : ncore(0, 1) = -1
4224 : END SELECT
4225 : END IF
4226 : END IF
4227 : ! if the core is established, finish the setup
4228 31291 : IF (ncore(0, 1) >= 0) THEN
4229 156455 : DO l = 0, lmin
4230 125164 : ll = 2*(2*l + 1)
4231 1376804 : nn = SUM(ncore(l, :)) + econfx(l)
4232 125164 : ii = 0
4233 31291 : DO
4234 149972 : ii = ii + 1
4235 149972 : IF (nn <= ll) THEN
4236 125164 : nelem(l, ii) = nn
4237 : EXIT
4238 : ELSE
4239 24808 : nelem(l, ii) = ll
4240 24808 : nn = nn - ll
4241 : END IF
4242 : END DO
4243 : END DO
4244 2221661 : ncalc = nelem - ncore
4245 : ELSE
4246 : ! test for compatibility of valence occupation and full atomic occupation
4247 0 : IF (iounit > 0) THEN
4248 0 : WRITE (iounit, "(/,A,A2)") "WARNING: Core states irregular for atom type ", ptable(z)%symbol
4249 0 : WRITE (iounit, "(A,10I3)") "WARNING: Redefine ELEC_CONF in the KIND section"
4250 0 : CPABORT("Incompatible Atomic Occupations Detected")
4251 : END IF
4252 : END IF
4253 : ELSE
4254 34 : lmin = MIN(lmat, UBOUND(ptable(z)%e_conv, 1))
4255 34 : ncore = 0
4256 34 : ncalc = 0
4257 170 : DO l = 0, lmin
4258 136 : ll = 2*(2*l + 1)
4259 136 : nn = ABS(econfx(l))
4260 136 : ii = 0
4261 34 : DO
4262 136 : ii = ii + 1
4263 136 : IF (nn <= ll) THEN
4264 136 : ncalc(l, ii) = -nn
4265 : EXIT
4266 : ELSE
4267 0 : ncalc(l, ii) = -ll
4268 0 : nn = nn - ll
4269 : END IF
4270 : END DO
4271 : END DO
4272 34 : nelem = ncalc
4273 : END IF
4274 :
4275 31325 : END SUBROUTINE set_pseudo_state
4276 :
4277 : ! **************************************************************************************************
4278 : !> \brief finds if a given qs run needs to use nlcc
4279 : !> \param qs_kind_set ...
4280 : !> \return ...
4281 : ! **************************************************************************************************
4282 55730 : FUNCTION has_nlcc(qs_kind_set) RESULT(nlcc)
4283 :
4284 : TYPE(qs_kind_type), DIMENSION(:) :: qs_kind_set
4285 : LOGICAL :: nlcc
4286 :
4287 : INTEGER :: ikind
4288 : LOGICAL :: nlcc_present
4289 : TYPE(gth_potential_type), POINTER :: gth_potential
4290 : TYPE(sgp_potential_type), POINTER :: sgp_potential
4291 :
4292 55730 : nlcc = .FALSE.
4293 :
4294 161731 : DO ikind = 1, SIZE(qs_kind_set)
4295 106001 : CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential, sgp_potential=sgp_potential)
4296 161731 : IF (ASSOCIATED(gth_potential)) THEN
4297 67755 : CALL get_potential(potential=gth_potential, nlcc_present=nlcc_present)
4298 67755 : nlcc = nlcc .OR. nlcc_present
4299 38246 : ELSE IF (ASSOCIATED(sgp_potential)) THEN
4300 1046 : CALL get_potential(potential=sgp_potential, has_nlcc=nlcc_present)
4301 1046 : nlcc = nlcc .OR. nlcc_present
4302 : END IF
4303 : END DO
4304 :
4305 55730 : END FUNCTION has_nlcc
4306 :
4307 : ! **************************************************************************************************
4308 :
4309 0 : END MODULE qs_kind_types
|
