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