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