(git:ccc2433)
qs_kind_types.F
Go to the documentation of this file.
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 ! **************************************************************************************************
27  USE atom_types, ONLY: atom_ecppot_type,&
28  lmat,&
30  USE atom_upf, ONLY: atom_read_upf,&
33  USE atomic_kind_types, ONLY: atomic_kind_type,&
36  basis_set_container_type,&
40  USE basis_set_types, ONLY: &
43  gto_basis_set_type, init_aux_basis_set, init_orb_basis_set, read_gto_basis_set, &
45  USE cp_control_types, ONLY: dft_control_type,&
46  qs_control_type,&
47  xtb_control_type
50  cp_logger_type
51  USE cp_output_handling, ONLY: cp_p_file,&
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
60  USE input_constants, ONLY: &
67  section_vals_type,&
69  USE kinds, ONLY: default_path_length,&
71  dp
72  USE mathconstants, ONLY: pi
73  USE message_passing, ONLY: mp_para_env_type
75  nco,&
76  ncoset
80  paw_proj_set_type,&
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
91  USE qs_dispersion_types, ONLY: qs_atom_dispersion_type
92  USE qs_grid_atom, ONLY: allocate_grid_atom,&
94  grid_atom_type
97  harmonics_atom_type
101  semi_empirical_type,&
106  USE string_utilities, ONLY: uppercase
107  USE xtb_parameters, ONLY: xtb_set_kab
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, &
238  get_qs_kind, &
239  get_qs_kind_set, &
240  has_nlcc, &
243  init_gapw_nlcc, &
245  set_qs_kind, &
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  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  IF (ASSOCIATED(qs_kind_set)) THEN
269 
270  nkind = SIZE(qs_kind_set)
271 
272  DO ikind = 1, nkind
273  IF (ASSOCIATED(qs_kind_set(ikind)%all_potential)) THEN
274  CALL deallocate_potential(qs_kind_set(ikind)%all_potential)
275  END IF
276  IF (ASSOCIATED(qs_kind_set(ikind)%tnadd_potential)) THEN
277  CALL deallocate_potential(qs_kind_set(ikind)%tnadd_potential)
278  END IF
279  IF (ASSOCIATED(qs_kind_set(ikind)%gth_potential)) THEN
280  CALL deallocate_potential(qs_kind_set(ikind)%gth_potential)
281  END IF
282  IF (ASSOCIATED(qs_kind_set(ikind)%sgp_potential)) THEN
283  CALL deallocate_potential(qs_kind_set(ikind)%sgp_potential)
284  END IF
285  IF (ASSOCIATED(qs_kind_set(ikind)%upf_potential)) THEN
286  CALL atom_release_upf(qs_kind_set(ikind)%upf_potential)
287  DEALLOCATE (qs_kind_set(ikind)%upf_potential)
288  END IF
289  IF (ASSOCIATED(qs_kind_set(ikind)%se_parameter)) THEN
290  CALL semi_empirical_release(qs_kind_set(ikind)%se_parameter)
291  END IF
292  IF (ASSOCIATED(qs_kind_set(ikind)%dftb_parameter)) THEN
293  CALL deallocate_dftb_atom_param(qs_kind_set(ikind)%dftb_parameter)
294  END IF
295  IF (ASSOCIATED(qs_kind_set(ikind)%xtb_parameter)) THEN
296  CALL deallocate_xtb_atom_param(qs_kind_set(ikind)%xtb_parameter)
297  END IF
298  IF (ASSOCIATED(qs_kind_set(ikind)%paw_proj_set)) THEN
299  CALL deallocate_paw_proj_set(qs_kind_set(ikind)%paw_proj_set)
300  END IF
301  IF (ASSOCIATED(qs_kind_set(ikind)%harmonics)) THEN
302  CALL deallocate_harmonics_atom(qs_kind_set(ikind)%harmonics)
303  END IF
304  IF (ASSOCIATED(qs_kind_set(ikind)%grid_atom)) THEN
305  CALL deallocate_grid_atom(qs_kind_set(ikind)%grid_atom)
306  END IF
307  IF (ASSOCIATED(qs_kind_set(ikind)%elec_conf)) THEN
308  DEALLOCATE (qs_kind_set(ikind)%elec_conf)
309  END IF
310 
311  IF (ASSOCIATED(qs_kind_set(ikind)%dft_plus_u)) THEN
312  IF (ASSOCIATED(qs_kind_set(ikind)%dft_plus_u%orbitals)) THEN
313  DEALLOCATE (qs_kind_set(ikind)%dft_plus_u%orbitals)
314  END IF
315  IF (ASSOCIATED(qs_kind_set(ikind)%dft_plus_u%nelec)) THEN
316  DEALLOCATE (qs_kind_set(ikind)%dft_plus_u%nelec)
317  END IF
318  DEALLOCATE (qs_kind_set(ikind)%dft_plus_u)
319  END IF
320 
321  IF (ASSOCIATED(qs_kind_set(ikind)%nlcc_pot)) THEN
322  DEALLOCATE (qs_kind_set(ikind)%nlcc_pot)
323  END IF
324 
325  IF (ASSOCIATED(qs_kind_set(ikind)%dispersion)) THEN
326  DEALLOCATE (qs_kind_set(ikind)%dispersion)
327  END IF
328  IF (ASSOCIATED(qs_kind_set(ikind)%addel)) THEN
329  DEALLOCATE (qs_kind_set(ikind)%addel)
330  END IF
331  IF (ASSOCIATED(qs_kind_set(ikind)%naddel)) THEN
332  DEALLOCATE (qs_kind_set(ikind)%naddel)
333  END IF
334  IF (ASSOCIATED(qs_kind_set(ikind)%laddel)) THEN
335  DEALLOCATE (qs_kind_set(ikind)%laddel)
336  END IF
337  IF (ASSOCIATED(qs_kind_set(ikind)%reltmat)) THEN
338  DEALLOCATE (qs_kind_set(ikind)%reltmat)
339  END IF
340 
341  IF (ASSOCIATED(qs_kind_set(ikind)%pao_potentials)) THEN
342  DEALLOCATE (qs_kind_set(ikind)%pao_potentials)
343  END IF
344  IF (ASSOCIATED(qs_kind_set(ikind)%pao_descriptors)) THEN
345  DEALLOCATE (qs_kind_set(ikind)%pao_descriptors)
346  END IF
347 
348  CALL remove_basis_set_container(qs_kind_set(ikind)%basis_sets)
349 
350  END DO
351  DEALLOCATE (qs_kind_set)
352  ELSE
353  CALL cp_abort(__location__, &
354  "The pointer qs_kind_set is not associated and "// &
355  "cannot be deallocated")
356  END IF
357 
358  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  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  IF (PRESENT(basis_type)) THEN
513  my_basis_type = basis_type
514  ELSE
515  my_basis_type = "ORB"
516  END IF
517 
518  IF (PRESENT(basis_set)) THEN
519  CALL get_basis_from_container(qs_kind%basis_sets, basis_set=basis_set, &
520  basis_type=my_basis_type)
521  END IF
522 
523  IF (PRESENT(ncgf)) THEN
524  CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis_set, &
525  basis_type=my_basis_type)
526  IF (ASSOCIATED(tmp_basis_set)) THEN
527  CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, ncgf=ncgf)
528  ELSE IF (ASSOCIATED(qs_kind%dftb_parameter)) THEN
529  l = qs_kind%dftb_parameter%lmax
530  ncgf = ((l + 1)*(l + 2)*(l + 3))/6
531  ELSE
532  ncgf = 0
533  END IF
534  END IF
535 
536  IF (PRESENT(nsgf)) THEN
537  CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis_set, &
538  basis_type=my_basis_type)
539  IF (ASSOCIATED(tmp_basis_set)) THEN
540  CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, nsgf=nsgf)
541  ELSE IF (ASSOCIATED(qs_kind%dftb_parameter)) THEN
542  nsgf = qs_kind%dftb_parameter%natorb
543  ELSE
544  nsgf = 0
545  END IF
546  END IF
547 
548  IF (PRESENT(all_potential)) all_potential => qs_kind%all_potential
549  IF (PRESENT(tnadd_potential)) tnadd_potential => qs_kind%tnadd_potential
550  IF (PRESENT(gth_potential)) gth_potential => qs_kind%gth_potential
551  IF (PRESENT(sgp_potential)) sgp_potential => qs_kind%sgp_potential
552  IF (PRESENT(upf_potential)) upf_potential => qs_kind%upf_potential
553  IF (PRESENT(se_parameter)) se_parameter => qs_kind%se_parameter
554  IF (PRESENT(dftb_parameter)) dftb_parameter => qs_kind%dftb_parameter
555  IF (PRESENT(xtb_parameter)) xtb_parameter => qs_kind%xtb_parameter
556 
557  IF (PRESENT(element_symbol)) element_symbol = qs_kind%element_symbol
558  IF (PRESENT(name)) name = qs_kind%name
559  IF (PRESENT(dftb3_param)) dftb3_param = qs_kind%dudq_dftb3
560  IF (PRESENT(elec_conf)) elec_conf => qs_kind%elec_conf
561  IF (PRESENT(alpha_core_charge)) THEN
562  IF (ASSOCIATED(qs_kind%all_potential)) THEN
563  CALL get_potential(potential=qs_kind%all_potential, &
564  alpha_core_charge=alpha_core_charge)
565  ELSE IF (ASSOCIATED(qs_kind%gth_potential)) THEN
566  CALL get_potential(potential=qs_kind%gth_potential, &
567  alpha_core_charge=alpha_core_charge)
568  ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
569  CALL get_potential(potential=qs_kind%sgp_potential, &
570  alpha_core_charge=alpha_core_charge)
571  ELSE
572  alpha_core_charge = 1.0_dp
573  END IF
574  END IF
575  IF (PRESENT(ccore_charge)) THEN
576  IF (ASSOCIATED(qs_kind%all_potential)) THEN
577  CALL get_potential(potential=qs_kind%all_potential, &
578  ccore_charge=ccore_charge)
579  ELSE IF (ASSOCIATED(qs_kind%gth_potential)) THEN
580  CALL get_potential(potential=qs_kind%gth_potential, &
581  ccore_charge=ccore_charge)
582  ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
583  CALL get_potential(potential=qs_kind%sgp_potential, &
584  ccore_charge=ccore_charge)
585  ELSE IF (ASSOCIATED(qs_kind%upf_potential)) THEN
586  cpabort("UPF CCORE CHARGE RADIUS NOT AVAILABLE")
587  ELSE
588  ccore_charge = 0.0_dp
589  END IF
590  END IF
591  IF (PRESENT(core_charge_radius)) THEN
592  IF (ASSOCIATED(qs_kind%all_potential)) THEN
593  CALL get_potential(potential=qs_kind%all_potential, &
594  core_charge_radius=core_charge_radius)
595  ELSE IF (ASSOCIATED(qs_kind%gth_potential)) THEN
596  CALL get_potential(potential=qs_kind%gth_potential, &
597  core_charge_radius=core_charge_radius)
598  ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
599  CALL get_potential(potential=qs_kind%sgp_potential, &
600  core_charge_radius=core_charge_radius)
601  ELSE IF (ASSOCIATED(qs_kind%upf_potential)) THEN
602  cpabort("UPF CORE CHARGE RADIUS NOT AVAILABLE")
603  ELSE
604  core_charge_radius = 0.0_dp
605  END IF
606  END IF
607  IF (PRESENT(core_charge)) THEN
608  IF (ASSOCIATED(qs_kind%all_potential)) THEN
609  CALL get_potential(potential=qs_kind%all_potential, &
610  zeff=core_charge)
611  ELSE IF (ASSOCIATED(qs_kind%gth_potential)) THEN
612  CALL get_potential(potential=qs_kind%gth_potential, &
613  zeff=core_charge)
614  ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
615  CALL get_potential(potential=qs_kind%sgp_potential, &
616  zeff=core_charge)
617  ELSE IF (ASSOCIATED(qs_kind%upf_potential)) THEN
618  cpabort("UPF CORE CHARGE NOT AVAILABLE")
619  ELSE
620  core_charge = 0.0_dp
621  END IF
622  END IF
623 
624  IF (PRESENT(zeff)) THEN
625  IF (ASSOCIATED(qs_kind%all_potential)) THEN
626  CALL get_potential(potential=qs_kind%all_potential, zeff=zeff)
627  ELSE IF (ASSOCIATED(qs_kind%gth_potential)) THEN
628  CALL get_potential(potential=qs_kind%gth_potential, zeff=zeff)
629  ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
630  CALL get_potential(potential=qs_kind%sgp_potential, zeff=zeff)
631  ELSE IF (ASSOCIATED(qs_kind%upf_potential)) THEN
632  zeff = qs_kind%upf_potential%zion
633  ELSE
634  zeff = 0.0_dp
635  END IF
636  END IF
637 
638  IF (PRESENT(covalent_radius)) covalent_radius = qs_kind%covalent_radius
639  IF (PRESENT(vdw_radius)) vdw_radius = qs_kind%vdw_radius
640 
641  IF (PRESENT(paw_proj_set)) paw_proj_set => qs_kind%paw_proj_set
642  IF (PRESENT(paw_atom)) paw_atom = qs_kind%paw_atom
643  IF (PRESENT(gpw_r3d_rs_type_forced)) gpw_r3d_rs_type_forced = qs_kind%gpw_r3d_rs_type_forced
644  IF (PRESENT(hard_radius)) hard_radius = qs_kind%hard_radius
645  IF (PRESENT(hard0_radius)) hard0_radius = qs_kind%hard0_radius
646  IF (PRESENT(max_rad_local)) max_rad_local = qs_kind%max_rad_local
647  IF (PRESENT(harmonics)) harmonics => qs_kind%harmonics
648  IF (PRESENT(max_s_harm)) THEN
649  IF (ASSOCIATED(qs_kind%harmonics)) THEN
650  max_s_harm = qs_kind%harmonics%max_s_harm
651  ELSE
652  max_s_harm = 0
653  END IF
654  END IF
655  IF (PRESENT(max_iso_not0)) THEN
656  IF (ASSOCIATED(qs_kind%harmonics)) THEN
657  max_iso_not0 = qs_kind%harmonics%max_iso_not0
658  ELSE
659  max_iso_not0 = 0
660  END IF
661  END IF
662  IF (PRESENT(grid_atom)) grid_atom => qs_kind%grid_atom
663  IF (PRESENT(ngrid_ang)) ngrid_ang = qs_kind%ngrid_ang
664  IF (PRESENT(ngrid_rad)) ngrid_rad = qs_kind%ngrid_rad
665  IF (PRESENT(lmax_rho0)) lmax_rho0 = qs_kind%lmax_rho0
666  IF (PRESENT(ghost)) ghost = qs_kind%ghost
667  IF (PRESENT(floating)) floating = qs_kind%floating
668  IF (PRESENT(dft_plus_u_atom)) dft_plus_u_atom = ASSOCIATED(qs_kind%dft_plus_u)
669  IF (PRESENT(l_of_dft_plus_u)) THEN
670  IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
671  l_of_dft_plus_u = qs_kind%dft_plus_u%l
672  ELSE
673  l_of_dft_plus_u = -1
674  END IF
675  END IF
676  IF (PRESENT(n_of_dft_plus_u)) THEN
677  IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
678  n_of_dft_plus_u = qs_kind%dft_plus_u%n
679  ELSE
680  n_of_dft_plus_u = -1
681  END IF
682  END IF
683  IF (PRESENT(u_minus_j)) THEN
684  IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
685  u_minus_j = qs_kind%dft_plus_u%u_minus_j
686  ELSE
687  u_minus_j = 0.0_dp
688  END IF
689  END IF
690  IF (PRESENT(u_minus_j_target)) THEN
691  IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
692  u_minus_j_target = qs_kind%dft_plus_u%u_minus_j_target
693  ELSE
694  u_minus_j_target = 0.0_dp
695  END IF
696  END IF
697  IF (PRESENT(u_of_dft_plus_u)) THEN
698  IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
699  u_of_dft_plus_u = qs_kind%dft_plus_u%U
700  ELSE
701  u_of_dft_plus_u = 0.0_dp
702  END IF
703  END IF
704  IF (PRESENT(j_of_dft_plus_u)) THEN
705  IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
706  j_of_dft_plus_u = qs_kind%dft_plus_u%J
707  ELSE
708  j_of_dft_plus_u = 0.0_dp
709  END IF
710  END IF
711  IF (PRESENT(alpha_of_dft_plus_u)) THEN
712  IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
713  alpha_of_dft_plus_u = qs_kind%dft_plus_u%alpha
714  ELSE
715  alpha_of_dft_plus_u = 0.0_dp
716  END IF
717  END IF
718  IF (PRESENT(beta_of_dft_plus_u)) THEN
719  IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
720  beta_of_dft_plus_u = qs_kind%dft_plus_u%beta
721  ELSE
722  beta_of_dft_plus_u = 0.0_dp
723  END IF
724  END IF
725  IF (PRESENT(j0_of_dft_plus_u)) THEN
726  IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
727  j0_of_dft_plus_u = qs_kind%dft_plus_u%J0
728  ELSE
729  j0_of_dft_plus_u = 0.0_dp
730  END IF
731  END IF
732  IF (PRESENT(occupation_of_dft_plus_u)) THEN
733  IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
734  occupation_of_dft_plus_u = qs_kind%dft_plus_u%occupation
735  ELSE
736  occupation_of_dft_plus_u = -1.0_dp
737  END IF
738  END IF
739 
740  IF (PRESENT(init_u_ramping_each_scf)) THEN
741  IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
742  init_u_ramping_each_scf = qs_kind%dft_plus_u%init_u_ramping_each_scf
743  ELSE
744  init_u_ramping_each_scf = .false.
745  END IF
746  END IF
747  IF (PRESENT(u_ramping)) THEN
748  IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
749  u_ramping = qs_kind%dft_plus_u%u_ramping
750  ELSE
751  u_ramping = 0.0_dp
752  END IF
753  END IF
754  IF (PRESENT(eps_u_ramping)) THEN
755  IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
756  eps_u_ramping = qs_kind%dft_plus_u%eps_u_ramping
757  ELSE
758  eps_u_ramping = 1.0e-5_dp
759  END IF
760  END IF
761  IF (PRESENT(nelec)) THEN
762  NULLIFY (nelec)
763  IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
764  IF (ASSOCIATED(qs_kind%dft_plus_u%nelec)) THEN
765  nelec => qs_kind%dft_plus_u%nelec
766  END IF
767  END IF
768  END IF
769  IF (PRESENT(orbitals)) THEN
770  NULLIFY (orbitals)
771  IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
772  IF (ASSOCIATED(qs_kind%dft_plus_u%orbitals)) THEN
773  orbitals => qs_kind%dft_plus_u%orbitals
774  END IF
775  END IF
776  END IF
777  IF (PRESENT(eps_scf)) THEN
778  IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
779  eps_scf = qs_kind%dft_plus_u%eps_scf
780  ELSE
781  eps_scf = 1.0e30_dp
782  END IF
783  END IF
784  IF (PRESENT(max_scf)) THEN
785  IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
786  max_scf = qs_kind%dft_plus_u%max_scf
787  ELSE
788  max_scf = -1
789  END IF
790  END IF
791  IF (PRESENT(smear)) THEN
792  IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
793  smear = qs_kind%dft_plus_u%smear
794  ELSE
795  smear = .false.
796  END IF
797  END IF
798  IF (PRESENT(dispersion)) dispersion => qs_kind%dispersion
799  IF (PRESENT(bs_occupation)) bs_occupation = qs_kind%bs_occupation
800  IF (PRESENT(addel)) addel => qs_kind%addel
801  IF (PRESENT(laddel)) laddel => qs_kind%laddel
802  IF (PRESENT(naddel)) naddel => qs_kind%naddel
803 
804  IF (PRESENT(magnetization)) magnetization = qs_kind%magnetization
805 
806  IF (PRESENT(no_optimize)) no_optimize = qs_kind%no_optimize
807 
808  IF (PRESENT(reltmat)) reltmat => qs_kind%reltmat
809 
810  IF (PRESENT(mao)) mao = qs_kind%mao
811 
812  IF (PRESENT(lmax_dftb)) lmax_dftb = qs_kind%lmax_dftb
813 
814  IF (PRESENT(pao_basis_size)) pao_basis_size = qs_kind%pao_basis_size
815  IF (PRESENT(pao_potentials)) pao_potentials => qs_kind%pao_potentials
816  IF (PRESENT(pao_descriptors)) pao_descriptors => qs_kind%pao_descriptors
817  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  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  IF (PRESENT(basis_type)) THEN
895  my_basis_type = basis_type
896  ELSE
897  my_basis_type = "ORB"
898  END IF
899 
900  IF (ASSOCIATED(qs_kind_set)) THEN
901 
902  IF (PRESENT(maxcgf)) maxcgf = 0
903  IF (PRESENT(maxco)) maxco = 0
904  IF (PRESENT(maxco_proj)) maxco_proj = 0
905  IF (PRESENT(maxg_iso_not0)) maxg_iso_not0 = 0
906  IF (PRESENT(maxgtops)) maxgtops = 0
907  IF (PRESENT(maxlgto)) maxlgto = -1
908  IF (PRESENT(maxlppl)) maxlppl = -1
909  IF (PRESENT(maxlppnl)) maxlppnl = -1
910  IF (PRESENT(maxpol)) maxpol = -1
911  IF (PRESENT(maxlprj)) maxlprj = -1
912  IF (PRESENT(maxnset)) maxnset = 0
913  IF (PRESENT(maxppnl)) maxppnl = 0
914  IF (PRESENT(maxsgf)) maxsgf = 0
915  IF (PRESENT(maxsgf_set)) maxsgf_set = 0
916  IF (PRESENT(ncgf)) ncgf = 0
917  IF (PRESENT(nelectron)) nelectron = 0
918  IF (PRESENT(npgf)) npgf = 0
919  IF (PRESENT(nset)) nset = 0
920  IF (PRESENT(nsgf)) nsgf = 0
921  IF (PRESENT(nshell)) nshell = 0
922  IF (PRESENT(all_potential_present)) all_potential_present = .false.
923  IF (PRESENT(tnadd_potential_present)) tnadd_potential_present = .false.
924  IF (PRESENT(gth_potential_present)) gth_potential_present = .false.
925  IF (PRESENT(sgp_potential_present)) sgp_potential_present = .false.
926  IF (PRESENT(paw_atom_present)) paw_atom_present = .false.
927  IF (PRESENT(max_ngrid_rad)) max_ngrid_rad = 0
928  IF (PRESENT(max_sph_harm)) max_sph_harm = 0
929  IF (PRESENT(lmax_rho0)) lmax_rho0 = 0
930  IF (PRESENT(basis_rcut)) basis_rcut = 0.0_dp
931  IF (PRESENT(total_zeff_corr)) total_zeff_corr = 0.0_dp
932 
933  nkind = SIZE(qs_kind_set)
934  DO ikind = 1, nkind
935  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  lmax_rho0=lmax_rho0_kind)
949 
950  IF (PRESENT(maxlppl) .AND. ASSOCIATED(gth_potential)) THEN
951  CALL get_potential(potential=gth_potential, nexp_ppl=n)
952  maxlppl = max(maxlppl, 2*(n - 1))
953  ELSEIF (PRESENT(maxlppl) .AND. ASSOCIATED(sgp_potential)) THEN
954  CALL get_potential(potential=sgp_potential, nrloc=nrloc, ecp_semi_local=ecp_semi_local)
955  n = maxval(nrloc) - 2
956  maxlppl = max(maxlppl, 2*(n - 1))
957  IF (ecp_semi_local) THEN
958  CALL get_potential(potential=sgp_potential, sl_lmax=imax, nrpot=nrpot)
959  n = maxval(nrpot) - 2
960  n = 2*(n - 1) + imax
961  maxlppl = max(maxlppl, n)
962  END IF
963  END IF
964 
965  IF (PRESENT(maxlppnl) .AND. ASSOCIATED(gth_potential)) THEN
966  CALL get_potential(potential=gth_potential, lprj_ppnl_max=imax)
967  maxlppnl = max(maxlppnl, imax)
968  ELSEIF (PRESENT(maxlppnl) .AND. ASSOCIATED(sgp_potential)) THEN
969  CALL get_potential(potential=sgp_potential, lmax=imax)
970  maxlppnl = max(maxlppnl, imax)
971  END IF
972 
973  IF (PRESENT(maxpol) .AND. ASSOCIATED(tnadd_potential)) THEN
974  CALL get_potential(potential=tnadd_potential, npol=n)
975  maxpol = max(maxpol, 2*(n - 1))
976  END IF
977 
978  IF (PRESENT(maxco_proj) .AND. ASSOCIATED(paw_proj_set)) THEN
979  CALL get_paw_proj_set(paw_proj_set=paw_proj_set, ncgauprj=imax)
980  maxco_proj = max(maxco_proj, imax)
981  END IF
982 
983  IF (PRESENT(maxlprj) .AND. ASSOCIATED(paw_proj_set)) THEN
984  CALL get_paw_proj_set(paw_proj_set=paw_proj_set, maxl=imax)
985  maxlprj = max(maxlprj, imax)
986  END IF
987 
988  IF (PRESENT(maxppnl) .AND. ASSOCIATED(gth_potential)) THEN
989  CALL get_potential(potential=gth_potential, nppnl=imax)
990  maxppnl = max(maxppnl, imax)
991  ELSEIF (PRESENT(maxppnl) .AND. ASSOCIATED(sgp_potential)) THEN
992  CALL get_potential(potential=sgp_potential, nppnl=imax)
993  maxppnl = max(maxppnl, imax)
994  END IF
995 
996  CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis_set, &
997  basis_type=my_basis_type)
998 
999  IF (PRESENT(maxcgf)) THEN
1000  IF (ASSOCIATED(tmp_basis_set)) THEN
1001  CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, ncgf=imax)
1002  maxcgf = max(maxcgf, imax)
1003  ELSE IF (ASSOCIATED(qs_kind%dftb_parameter)) THEN
1004  CALL get_dftb_atom_param(dftb_parameter=dftb_parameter, lmax=imax)
1005  imax = ((imax + 1)*(imax + 2)*(imax + 3))/6
1006  maxcgf = max(maxcgf, imax)
1007  END IF
1008  END IF
1009 
1010  IF (PRESENT(maxco)) THEN
1011  IF (ASSOCIATED(tmp_basis_set)) THEN
1012  IF (PRESENT(maxder)) THEN
1013  CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, &
1014  maxco=imax, maxder=maxder)
1015  ELSE
1016  CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, maxco=imax)
1017  END IF
1018  maxco = max(maxco, imax)
1019  END IF
1020  IF (ASSOCIATED(gth_potential)) THEN
1021  CALL get_potential(potential=gth_potential, lprj_ppnl_max=imax)
1022  maxco = max(maxco, ncoset(imax))
1023  END IF
1024  IF (ASSOCIATED(sgp_potential)) THEN
1025  CALL get_potential(potential=sgp_potential, lmax=imax)
1026  maxco = max(maxco, ncoset(imax))
1027  CALL get_potential(potential=sgp_potential, sl_lmax=imax)
1028  maxco = max(maxco, ncoset(imax))
1029  END IF
1030  END IF
1031 
1032  IF (PRESENT(maxgtops)) THEN
1033  IF (ASSOCIATED(tmp_basis_set)) THEN
1034  CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, maxso=imax, nset=n)
1035  maxgtops = max(maxgtops, n*imax)
1036  END IF
1037  END IF
1038 
1039  IF (PRESENT(maxlgto)) THEN
1040  IF (ASSOCIATED(tmp_basis_set)) THEN
1041  CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, maxl=imax)
1042  maxlgto = max(maxlgto, imax)
1043  ELSE IF (ASSOCIATED(qs_kind%dftb_parameter)) THEN
1044  CALL get_dftb_atom_param(dftb_parameter=dftb_parameter, lmax=imax)
1045  maxlgto = max(maxlgto, imax)
1046  END IF
1047  END IF
1048 
1049  IF (PRESENT(maxnset)) THEN
1050  IF (ASSOCIATED(tmp_basis_set)) THEN
1051  CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, nset=n)
1052  maxnset = max(maxnset, n)
1053  END IF
1054  END IF
1055 
1056  IF (PRESENT(maxsgf)) THEN
1057  IF (ASSOCIATED(tmp_basis_set)) THEN
1058  CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, nsgf=imax)
1059  maxsgf = max(maxsgf, imax)
1060  END IF
1061  END IF
1062 
1063  IF (PRESENT(maxsgf_set)) THEN
1064  IF (ASSOCIATED(tmp_basis_set)) THEN
1065  CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, maxsgf_set=imax)
1066  maxsgf_set = max(maxsgf_set, imax)
1067  END IF
1068  END IF
1069 
1070  IF (PRESENT(ncgf)) THEN
1071  IF (ASSOCIATED(tmp_basis_set)) THEN
1072  CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, ncgf=n)
1073  ncgf = ncgf + n*qs_kind_set(ikind)%natom
1074  ELSE IF (ASSOCIATED(qs_kind%dftb_parameter)) THEN
1075  CALL get_dftb_atom_param(dftb_parameter=dftb_parameter, lmax=imax)
1076  n = ((imax + 1)*(imax + 2)*(imax + 3))/6
1077  ncgf = ncgf + n*qs_kind_set(ikind)%natom
1078  END IF
1079  END IF
1080 
1081  IF (PRESENT(npgf)) THEN
1082  IF (ASSOCIATED(tmp_basis_set)) THEN
1083  CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, npgf_sum=n)
1084  npgf = npgf + n*qs_kind_set(ikind)%natom
1085  END IF
1086  END IF
1087 
1088  IF (PRESENT(nset)) THEN
1089  IF (ASSOCIATED(tmp_basis_set)) THEN
1090  CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, nset=n)
1091  nset = nset + n*qs_kind_set(ikind)%natom
1092  END IF
1093  END IF
1094 
1095  IF (PRESENT(nsgf)) THEN
1096  IF (ASSOCIATED(tmp_basis_set)) THEN
1097  CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, nsgf=n)
1098  nsgf = nsgf + n*qs_kind_set(ikind)%natom
1099  ELSE IF (ASSOCIATED(qs_kind%dftb_parameter)) THEN
1100  CALL get_dftb_atom_param(dftb_parameter=dftb_parameter, natorb=n)
1101  nsgf = nsgf + n*qs_kind_set(ikind)%natom
1102  END IF
1103  END IF
1104 
1105  IF (PRESENT(nshell)) THEN
1106  IF (ASSOCIATED(tmp_basis_set)) THEN
1107  CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, nshell_sum=n)
1108  nshell = nshell + n*qs_kind_set(ikind)%natom
1109  ELSE IF (ASSOCIATED(qs_kind%dftb_parameter)) THEN
1110  CALL get_dftb_atom_param(dftb_parameter=dftb_parameter, lmax=n)
1111  nshell = nshell + (n + 1)*qs_kind_set(ikind)%natom
1112  END IF
1113  END IF
1114 
1115  IF (PRESENT(nelectron)) THEN
1116  IF (ASSOCIATED(qs_kind%all_potential)) THEN
1117  CALL get_potential(potential=qs_kind%all_potential, &
1118  zeff=zeff, zeff_correction=zeff_correction)
1119  ELSE IF (ASSOCIATED(qs_kind%gth_potential)) THEN
1120  CALL get_potential(potential=qs_kind%gth_potential, &
1121  zeff=zeff, zeff_correction=zeff_correction)
1122  ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
1123  CALL get_potential(potential=qs_kind%sgp_potential, &
1124  zeff=zeff, zeff_correction=zeff_correction)
1125  ELSE
1126  zeff = 0.0_dp
1127  zeff_correction = 0.0_dp
1128  END IF
1129  nelectron = nelectron + qs_kind_set(ikind)%natom*nint(zeff - zeff_correction)
1130  END IF
1131 
1132  IF (PRESENT(basis_rcut)) THEN
1133  IF (ASSOCIATED(tmp_basis_set)) THEN
1134  CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, kind_radius=brcut)
1135  basis_rcut = max(basis_rcut, brcut)
1136  ELSE IF (ASSOCIATED(qs_kind%dftb_parameter)) THEN
1137  CALL get_dftb_atom_param(dftb_parameter=dftb_parameter, cutoff=brcut)
1138  basis_rcut = max(basis_rcut, brcut)
1139  END IF
1140  END IF
1141 
1142  IF (PRESENT(total_zeff_corr)) THEN
1143  IF (ASSOCIATED(qs_kind%all_potential)) THEN
1144  CALL get_potential(potential=qs_kind%all_potential, &
1145  zeff=zeff, zeff_correction=zeff_correction)
1146  ELSE IF (ASSOCIATED(qs_kind%gth_potential)) THEN
1147  CALL get_potential(potential=qs_kind%gth_potential, &
1148  zeff=zeff, zeff_correction=zeff_correction)
1149  ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
1150  CALL get_potential(potential=qs_kind%sgp_potential, &
1151  zeff=zeff, zeff_correction=zeff_correction)
1152  ELSE
1153  zeff = 0.0_dp
1154  zeff_correction = 0.0_dp
1155  END IF
1156  total_zeff_corr = total_zeff_corr + qs_kind_set(ikind)%natom*zeff_correction
1157  END IF
1158 
1159  IF (PRESENT(all_potential_present)) THEN
1160  IF (ASSOCIATED(all_potential)) THEN
1161  all_potential_present = .true.
1162  END IF
1163  END IF
1164 
1165  IF (PRESENT(tnadd_potential_present)) THEN
1166  IF (ASSOCIATED(tnadd_potential)) THEN
1167  tnadd_potential_present = .true.
1168  END IF
1169  END IF
1170 
1171  IF (PRESENT(gth_potential_present)) THEN
1172  IF (ASSOCIATED(gth_potential)) THEN
1173  gth_potential_present = .true.
1174  END IF
1175  END IF
1176 
1177  IF (PRESENT(sgp_potential_present)) THEN
1178  IF (ASSOCIATED(sgp_potential)) THEN
1179  sgp_potential_present = .true.
1180  END IF
1181  END IF
1182 
1183  IF (PRESENT(paw_atom_present)) THEN
1184  IF (paw_atom) THEN
1185  paw_atom_present = .true.
1186  END IF
1187  END IF
1188 
1189  IF (PRESENT(dft_plus_u_atom_present)) THEN
1190  IF (dft_plus_u_atom) THEN
1191  dft_plus_u_atom_present = .true.
1192  END IF
1193  END IF
1194 
1195  IF (PRESENT(max_ngrid_rad)) THEN
1196  max_ngrid_rad = max(max_ngrid_rad, ngrid_rad)
1197  END IF
1198 
1199  IF (PRESENT(max_sph_harm)) THEN
1200  max_sph_harm = max(max_sph_harm, max_s_harm)
1201  END IF
1202 
1203  IF (PRESENT(maxg_iso_not0)) THEN
1204  maxg_iso_not0 = max(maxg_iso_not0, max_iso_not0)
1205  END IF
1206 
1207  IF (PRESENT(lmax_rho0)) THEN
1208  lmax_rho0 = max(lmax_rho0, lmax_rho0_kind)
1209  END IF
1210 
1211  END DO
1212  ELSE
1213  cpabort("The pointer qs_kind_set is not associated")
1214  END IF
1215 
1216  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  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  CALL timeset(routinen, handle)
1234 
1235  cpassert(ASSOCIATED(qs_kind))
1236 
1237  IF (ASSOCIATED(qs_kind%gth_potential)) THEN
1238  CALL init_potential(qs_kind%gth_potential)
1239  ELSEIF (ASSOCIATED(qs_kind%sgp_potential)) THEN
1240  CALL init_potential(qs_kind%sgp_potential)
1241  END IF
1242 
1243  DO i = 1, SIZE(qs_kind%basis_sets, 1)
1244  NULLIFY (tmp_basis_set)
1245  CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis_set, &
1246  inumbas=i, basis_type=basis_type)
1247  IF (basis_type == "") cycle
1248  IF (basis_type == "AUX") THEN
1249  IF (tmp_basis_set%norm_type < 0) tmp_basis_set%norm_type = 1
1250  CALL init_aux_basis_set(tmp_basis_set)
1251  ELSE
1252  IF (tmp_basis_set%norm_type < 0) tmp_basis_set%norm_type = 2
1253  CALL init_orb_basis_set(tmp_basis_set)
1254  END IF
1255  END DO
1256 
1257  CALL timestop(handle)
1258 
1259  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  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  CALL timeset(routinen, handle)
1277 
1278  IF (.NOT. ASSOCIATED(qs_kind_set)) THEN
1279  cpabort("init_qs_kind_set: The pointer qs_kind_set is not associated")
1280  END IF
1281 
1282  DO ikind = 1, SIZE(qs_kind_set)
1283  qs_kind => qs_kind_set(ikind)
1284  CALL init_qs_kind(qs_kind)
1285  END DO
1286 
1287  CALL timestop(handle)
1288 
1289  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  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  my_mod_control = .true.
1314  IF (PRESENT(modify_qs_control)) THEN
1315  my_mod_control = modify_qs_control
1316  END IF
1317 
1318  IF (ASSOCIATED(qs_kind_set)) THEN
1319 
1320  IF (my_mod_control) qs_control%gapw_control%non_paw_atoms = .false.
1321  nkind = SIZE(qs_kind_set)
1322 
1323  DO ikind = 1, nkind
1324 
1325  qs_kind => qs_kind_set(ikind)
1326 
1327  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  max_rad_local=max_rad_local_type, gpw_r3d_rs_type_forced=gpw)
1330 
1331  NULLIFY (soft_basis)
1332  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  qs_control%gapw_control%force_paw, gpw)
1336  CALL add_basis_set_to_container(qs_kind%basis_sets, soft_basis, "ORB_SOFT")
1337  CALL set_qs_kind(qs_kind=qs_kind, paw_atom=paw_atom)
1338 
1339  bas1c = qs_control%gapw_control%basis_1c
1340  NULLIFY (basis_1c)
1341  SELECT CASE (bas1c)
1342  CASE (gapw_1c_orb)
1343  ilevel = 0
1344  CASE (gapw_1c_small)
1345  ilevel = 1
1346  CASE (gapw_1c_medium)
1347  ilevel = 2
1348  CASE (gapw_1c_large)
1349  ilevel = 3
1350  CASE (gapw_1c_very_large)
1351  ilevel = 4
1352  CASE DEFAULT
1353  cpabort("basis_1c type")
1354  END SELECT
1355  CALL remove_basis_from_container(qs_kind%basis_sets, basis_type="GAPW_1C")
1356  CALL create_1c_basis(orb_basis, soft_basis, basis_1c, ilevel)
1357  CALL get_gto_basis_set(gto_basis_set=orb_basis, name=bsname)
1358  basis_1c%name = trim(bsname)//"_1c"
1359  CALL add_basis_set_to_container(qs_kind%basis_sets, basis_1c, "GAPW_1C")
1360  IF (paw_atom) THEN
1361  CALL allocate_paw_proj_set(qs_kind%paw_proj_set)
1362  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  max_rad_local_type, force_env_section)
1365  ELSE
1366  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  NULLIFY (qs_kind%grid_atom, qs_kind%harmonics)
1371  CALL allocate_grid_atom(qs_kind%grid_atom)
1372  CALL allocate_harmonics_atom(qs_kind%harmonics)
1373 
1374  END DO
1375 
1376  IF (my_mod_control) THEN
1377  IF (qs_control%gapw_control%non_paw_atoms) THEN
1378  qs_control%gapw_control%nopaw_as_gpw = .true.
1379  ELSE
1380  qs_control%gapw_control%nopaw_as_gpw = .false.
1381  END IF
1382  END IF
1383  ELSE
1384  cpabort("The pointer qs_kind_set is not associated")
1385  END IF
1386 
1387  END SUBROUTINE init_gapw_basis_set
1388 ! **************************************************************************************************
1389 !> \brief ...
1390 !> \param qs_kind_set ...
1391 ! **************************************************************************************************
1392  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  INTEGER, DIMENSION(:), POINTER :: nct_nlcc
1399  LOGICAL :: nlcc, nlcc_type, paw_atom
1400  REAL(dp) :: alpha, coa, cval
1401  REAL(kind=dp), DIMENSION(:), POINTER :: a_nlcc, alpha_nlcc, c_nlcc, fe, rc, rr
1402  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  IF (ASSOCIATED(qs_kind_set)) THEN
1408  nlcc = has_nlcc(qs_kind_set)
1409  IF (nlcc) THEN
1410  nkind = SIZE(qs_kind_set)
1411  DO ikind = 1, nkind
1412  qs_kind => qs_kind_set(ikind)
1413  CALL get_qs_kind(qs_kind, paw_atom=paw_atom)
1414  IF (paw_atom) THEN
1415  CALL get_qs_kind(qs_kind, gth_potential=gth_potential)
1416  CALL get_qs_kind(qs_kind, sgp_potential=sgp_potential)
1417  IF (ASSOCIATED(gth_potential)) THEN
1418  CALL get_potential(potential=gth_potential, nlcc_present=nlcc_type, &
1419  nexp_nlcc=nexp_nlcc, alpha_nlcc=alpha_nlcc, nct_nlcc=nct_nlcc, cval_nlcc=cval_nlcc)
1420  IF (nlcc_type) THEN
1421  nr = qs_kind%grid_atom%nr
1422  rr => qs_kind%grid_atom%rad
1423  ALLOCATE (qs_kind%nlcc_pot(nr, 2), rc(nr), fe(nr))
1424  den => qs_kind%nlcc_pot
1425  den = 0.0_dp
1426  DO i = 1, nexp_nlcc
1427  alpha = alpha_nlcc(i)
1428  rc(:) = rr(:)/alpha
1429  fe(:) = exp(-0.5_dp*rc(:)*rc(:))
1430  nc = nct_nlcc(i)
1431  DO ic = 1, nc
1432  cval = cval_nlcc(ic, i)
1433  coa = cval/alpha
1434  den(:, 1) = den(:, 1) + fe(:)*rc**(2*ic - 2)*cval
1435  den(:, 2) = den(:, 2) - fe(:)*rc**(2*ic - 1)*coa
1436  IF (ic > 1) THEN
1437  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  DEALLOCATE (rc, fe)
1442  END IF
1443  ELSE IF (ASSOCIATED(sgp_potential)) THEN
1444  CALL get_potential(potential=sgp_potential, has_nlcc=nlcc_type, &
1445  n_nlcc=n_nlcc, a_nlcc=a_nlcc, c_nlcc=c_nlcc)
1446  IF (nlcc_type) THEN
1447  nr = qs_kind%grid_atom%nr
1448  rr => qs_kind%grid_atom%rad
1449  ALLOCATE (qs_kind%nlcc_pot(nr, 2), rc(nr), fe(nr))
1450  den => qs_kind%nlcc_pot
1451  den = 0.0_dp
1452  DO i = 1, n_nlcc
1453  alpha = a_nlcc(i)
1454  fe(:) = exp(-alpha*rr(:)*rr(:))
1455  cval = c_nlcc(i)
1456  den(:, 1) = den(:, 1) + cval*fe(:)
1457  den(:, 2) = den(:, 2) - 2.0_dp*alpha*cval*rr(:)*fe(:)
1458  END DO
1459  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  cpabort("The pointer qs_kind_set is not associated")
1469  END IF
1470 
1471  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  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  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  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  REAL(kind=dp), DIMENSION(:), POINTER :: a_nl, aloc, anlcc, cloc, cnlcc, nelec
1519  REAL(kind=dp), DIMENSION(:, :), POINTER :: h_nl
1520  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: c_nl
1521  TYPE(atom_ecppot_type) :: ecppot
1522  TYPE(atom_sgp_potential_type) :: sgppot
1523  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  CALL timeset(routinen, handle)
1533 
1534  NULLIFY (logger)
1535  logger => cp_get_default_logger()
1536  iounit = cp_logger_get_default_io_unit(logger)
1537 
1538  NULLIFY (elec_conf)
1539 
1540  update_input = .true.
1541  basis_set_name(:) = ""
1542  basis_set_type(:) = ""
1543  basis_set_form(:) = ""
1544  potential_name = ""
1545  potential_type = ""
1546  kgpot_name = ""
1547  kgpot_type = ""
1548  z = -1
1549  zeff_correction = 0.0_dp
1550  explicit = .false.
1551  explicit_basis = .false.
1552  explicit_j = .false.
1553  explicit_kgpot = .false.
1554  explicit_potential = .false.
1555  explicit_u = .false.
1556  explicit_u_m_j = .false.
1557 
1558  dft_section => section_vals_get_subs_vals(force_env_section, "DFT")
1559  CALL section_vals_get(kind_section, n_repetition=n_rep)
1560  k_rep = -1
1561  akind_name = qs_kind%name
1562  CALL uppercase(akind_name)
1563  ! First we use the atom_name to find out the proper KIND section
1564  DO i_rep = 1, n_rep
1565  CALL section_vals_val_get(kind_section, "_SECTION_PARAMETERS_", &
1566  c_val=keyword, i_rep_section=i_rep)
1567  CALL uppercase(keyword)
1568  IF (keyword == akind_name) THEN
1569  k_rep = i_rep
1570  EXIT
1571  END IF
1572  END DO
1573  ! The search for the KIND section failed.. check for a QM/MM link atom
1574  IF (k_rep < 1) THEN
1575  ipos = index(qs_kind%name, "_")
1576  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  akind_name = qs_kind%name(1:ipos - 1)
1580  CALL uppercase(akind_name)
1581  DO i_rep = 1, n_rep
1582  CALL section_vals_val_get(kind_section, "_SECTION_PARAMETERS_", &
1583  c_val=keyword, i_rep_section=i_rep)
1584  CALL uppercase(keyword)
1585  IF (keyword == akind_name) THEN
1586  k_rep = i_rep
1587  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  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  element_symbol = qs_kind%element_symbol(1:2)
1597  CALL uppercase(element_symbol)
1598  DO i_rep = 1, n_rep
1599  CALL section_vals_val_get(kind_section, "_SECTION_PARAMETERS_", &
1600  c_val=keyword, i_rep_section=i_rep)
1601  CALL uppercase(keyword)
1602  IF (keyword == element_symbol) THEN
1603  k_rep = i_rep
1604  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  IF (k_rep < 1) THEN
1611  DO i_rep = 1, n_rep
1612  CALL section_vals_val_get(kind_section, "_SECTION_PARAMETERS_", &
1613  c_val=keyword, i_rep_section=i_rep)
1614  CALL uppercase(keyword)
1615  IF (keyword == "DEFAULT") THEN
1616  update_input = .false.
1617  k_rep = i_rep
1618  EXIT
1619  END IF
1620  END DO
1621  END IF
1622  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  " and for the DEFAULT section but no match was found. Check your input file!")
1628  END IF
1629  ! Retrieve information on element
1630  CALL get_ptable_info(qs_kind%element_symbol, ielement=z)
1631 
1632  ! Normal parsing of the KIND section
1633  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  n_rep_val=nb_rep)
1639  IF (.NOT. explicit) nb_rep = 0
1640  cpassert(nb_rep <= maxbas)
1641  DO i = 1, nb_rep
1642  CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1643  keyword_name="BASIS_SET", i_rep_val=i, c_vals=tmpstringlist)
1644  IF (SIZE(tmpstringlist) == 1) THEN
1645  ! default is orbital type and GTO
1646  basis_set_type(i) = "ORB"
1647  basis_set_form(i) = "GTO"
1648  basis_set_name(i) = tmpstringlist(1)
1649  ELSEIF (SIZE(tmpstringlist) == 2) THEN
1650  ! default is GTO
1651  basis_set_type(i) = tmpstringlist(1)
1652  basis_set_form(i) = "GTO"
1653  basis_set_name(i) = tmpstringlist(2)
1654  ELSEIF (SIZE(tmpstringlist) == 3) THEN
1655  basis_set_type(i) = tmpstringlist(1)
1656  basis_set_form(i) = tmpstringlist(2)
1657  basis_set_name(i) = tmpstringlist(3)
1658  ELSE
1659  CALL cp_abort(__location__, &
1660  "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  IF (basis_set_form(i) /= "GTO" .AND. basis_set_form(i) /= "STO") THEN
1664  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  i_val=qs_kind%pao_basis_size)
1671 
1672  ! parse PAO_POTENTIAL sections
1673  pao_pot_section => section_vals_get_subs_vals(kind_section, "PAO_POTENTIAL", i_rep_section=k_rep)
1674  CALL section_vals_get(pao_pot_section, n_repetition=npaopot)
1675  ALLOCATE (qs_kind%pao_potentials(npaopot))
1676  DO ipaopot = 1, npaopot
1677  CALL section_vals_val_get(pao_pot_section, keyword_name="MAXL", i_rep_section=ipaopot, &
1678  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  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  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  r_val=qs_kind%pao_potentials(ipaopot)%weight)
1685  END DO
1686 
1687  ! parse PAO_DESCRIPTOR sections
1688  pao_desc_section => section_vals_get_subs_vals(kind_section, "PAO_DESCRIPTOR", i_rep_section=k_rep)
1689  CALL section_vals_get(pao_desc_section, n_repetition=npaodesc)
1690  ALLOCATE (qs_kind%pao_descriptors(npaodesc))
1691  DO ipaodesc = 1, npaodesc
1692  CALL section_vals_val_get(pao_desc_section, keyword_name="BETA", i_rep_section=ipaodesc, &
1693  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  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  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  keyword_name="ELEC_CONF", n_rep_val=i)
1703  IF (i > 0) THEN
1704  CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1705  keyword_name="ELEC_CONF", i_vals=elec_conf)
1706  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  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  keyword_name="POTENTIAL_FILE_NAME", c_val=potential_fn_kind)
1713  CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1714  keyword_name="POTENTIAL_TYPE", c_val=potential_type)
1715  CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1716  explicit=explicit, keyword_name="POTENTIAL", c_vals=tmpstringlist)
1717  IF (explicit) THEN
1718  IF (SIZE(tmpstringlist) == 1) THEN
1719  ! old type of input: start of name defines type
1720  potential_name = tmpstringlist(1)
1721  IF (potential_type == "") THEN
1722  ipos = index(potential_name, "-")
1723  IF (ipos > 1) THEN
1724  potential_type = potential_name(:ipos - 1)
1725  ELSE
1726  potential_type = potential_name
1727  END IF
1728  END IF
1729  ELSEIF (SIZE(tmpstringlist) == 2) THEN
1730  potential_type = tmpstringlist(1)
1731  potential_name = tmpstringlist(2)
1732  ELSE
1733  cpabort("POTENTIAL input list is not correct")
1734  END IF
1735  END IF
1736  CALL uppercase(potential_type)
1737 
1738  ! Parse KG POTENTIAL
1739  CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1740  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  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  keyword_name="ECP_SEMI_LOCAL", l_val=ecp_semi_local)
1747 
1748  ! Assign atomic covalent radius
1749  qs_kind%covalent_radius = ptable(z)%covalent_radius*bohr
1750  CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1751  keyword_name="COVALENT_RADIUS", r_val=r)
1752  IF (r > 0.0_dp) qs_kind%covalent_radius = r
1753 
1754  ! Assign atomic van der Waals radius
1755  qs_kind%vdw_radius = ptable(z)%vdw_radius*bohr
1756  CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1757  keyword_name="VDW_RADIUS", r_val=r)
1758  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  keyword_name="HARD_EXP_RADIUS")
1763  IF (i == 0) THEN
1764  IF (z == 1) THEN
1765  qs_kind%hard_radius = 1.2_dp
1766  ELSE
1767  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  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  keyword_name="RHO0_EXP_RADIUS")
1777  IF (i == 0) THEN
1778  qs_kind%hard0_radius = qs_kind%hard_radius
1779  ELSE
1780  CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1781  keyword_name="RHO0_EXP_RADIUS", r_val=qs_kind%hard0_radius)
1782  END IF
1783  IF (qs_kind%hard_radius < qs_kind%hard0_radius) &
1784  cpabort("rc0 should be <= rc")
1785 
1786  CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1787  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  keyword_name="LEBEDEV_GRID", i_val=qs_kind%ngrid_ang)
1790  IF (qs_kind%ngrid_ang <= 0) &
1791  cpabort("# point lebedev grid < 0")
1792  CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1793  keyword_name="RADIAL_GRID", i_val=qs_kind%ngrid_rad)
1794  IF (qs_kind%ngrid_rad <= 0) &
1795  cpabort("# point radial grid < 0")
1796  CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1797  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  keyword_name="GHOST", l_val=qs_kind%ghost)
1800  CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1801  keyword_name="FLOATING_BASIS_CENTER", l_val=qs_kind%floating)
1802  CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1803  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  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  keyword_name="DFTB3_PARAM", r_val=qs_kind%dudq_dftb3)
1811  CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
1812  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  keyword_name="MAO", i_val=qs_kind%mao)
1817 
1818  ! Read the BS subsection of the current atomic kind, if enabled
1819  NULLIFY (bs_section)
1820  bs_section => section_vals_get_subs_vals(kind_section, "BS", &
1821  i_rep_section=k_rep)
1822  section_enabled = .false.
1823  CALL section_vals_val_get(bs_section, "_SECTION_PARAMETERS_", &
1824  l_val=section_enabled)
1825  IF (section_enabled) THEN
1826  ! test for conflict with magnetization
1827  IF (qs_kind%magnetization /= 0.0_dp) THEN
1828  CALL cp_abort(__location__, "BS Section is in conflict with non-zero magnetization "// &
1829  "for this atom kind.")
1830  END IF
1831  qs_kind%bs_occupation = .true.
1832  !Alpha spin
1833  NULLIFY (spin_section)
1834  spin_section => section_vals_get_subs_vals(bs_section, "ALPHA")
1835  CALL section_vals_get(spin_section, explicit=explicit)
1836  IF (explicit) THEN
1837  NULLIFY (add_el)
1838  CALL section_vals_val_get(spin_section, &
1839  keyword_name="NEL", i_vals=add_el)
1840  cpassert(ASSOCIATED(add_el))
1841  ALLOCATE (qs_kind%addel(SIZE(add_el), 2))
1842  qs_kind%addel = 0
1843  qs_kind%addel(1:SIZE(add_el), 1) = add_el(1:SIZE(add_el))
1844  NULLIFY (add_el)
1845  CALL section_vals_val_get(spin_section, &
1846  keyword_name="L", i_vals=add_el)
1847  cpassert(ASSOCIATED(add_el))
1848  cpassert(SIZE(add_el) == SIZE(qs_kind%addel, 1))
1849  ALLOCATE (qs_kind%laddel(SIZE(add_el), 2))
1850  qs_kind%laddel = 0
1851  qs_kind%laddel(1:SIZE(add_el), 1) = add_el(1:SIZE(add_el))
1852  ALLOCATE (qs_kind%naddel(SIZE(add_el), 2))
1853  qs_kind%naddel = 0
1854  NULLIFY (add_el)
1855  CALL section_vals_val_get(spin_section, &
1856  keyword_name="N", n_rep_val=i)
1857  IF (i > 0) THEN
1858  CALL section_vals_val_get(spin_section, &
1859  keyword_name="N", i_vals=add_el)
1860  IF (SIZE(add_el) == SIZE(qs_kind%addel, 1)) THEN
1861  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  NULLIFY (spin_section)
1867  spin_section => section_vals_get_subs_vals(bs_section, "BETA")
1868  CALL section_vals_get(spin_section, explicit=explicit)
1869  IF (explicit) THEN
1870  NULLIFY (add_el)
1871  CALL section_vals_val_get(spin_section, &
1872  keyword_name="NEL", i_vals=add_el)
1873  cpassert(SIZE(add_el) == SIZE(qs_kind%addel, 1))
1874  qs_kind%addel(1:SIZE(add_el), 2) = add_el(1:SIZE(add_el))
1875  qs_kind%addel(:, :) = qs_kind%addel(:, :)
1876  NULLIFY (add_el)
1877  CALL section_vals_val_get(spin_section, &
1878  keyword_name="L", i_vals=add_el)
1879  cpassert(SIZE(add_el) == SIZE(qs_kind%addel, 1))
1880  qs_kind%laddel(1:SIZE(add_el), 2) = add_el(1:SIZE(add_el))
1881 
1882  CALL section_vals_val_get(spin_section, &
1883  keyword_name="N", n_rep_val=i)
1884  IF (i > 0) THEN
1885  NULLIFY (add_el)
1886  CALL section_vals_val_get(spin_section, &
1887  keyword_name="N", i_vals=add_el)
1888  IF (SIZE(add_el) == SIZE(qs_kind%addel, 1)) THEN
1889  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  NULLIFY (dft_plus_u_section)
1898  dft_plus_u_section => section_vals_get_subs_vals(kind_section, &
1899  subsection_name="DFT_PLUS_U", &
1900  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  l_val=section_enabled)
1905  IF (section_enabled) THEN
1906  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  i_val=l)
1912  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  i_val=nu)
1917  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  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  explicit=explicit_j)
1928 
1929  CALL section_vals_val_get(dft_plus_u_section, &
1930  keyword_name="alpha", &
1931  r_val=qs_kind%dft_plus_u%alpha)
1932 
1933  CALL section_vals_val_get(dft_plus_u_section, &
1934  keyword_name="beta", &
1935  r_val=qs_kind%dft_plus_u%beta)
1936 
1937  CALL section_vals_val_get(dft_plus_u_section, &
1938  keyword_name="J0", &
1939  r_val=qs_kind%dft_plus_u%J0)
1940 
1941  CALL section_vals_val_get(dft_plus_u_section, &
1942  keyword_name="occupation", &
1943  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  explicit=explicit_u_m_j)
1952 
1953  IF ((explicit_u .OR. explicit_j) .AND. explicit_u_m_j) THEN
1954  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  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  l_val=qs_kind%dft_plus_u%init_u_ramping_each_scf)
1963  IF (qs_kind%dft_plus_u%u_ramping > 0.0_dp) THEN
1964  qs_kind%dft_plus_u%u_minus_j = 0.0_dp
1965  ELSE
1966  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  r_val=qs_kind%dft_plus_u%eps_u_ramping)
1971 
1972  NULLIFY (enforce_occupation_section)
1973  enforce_occupation_section => section_vals_get_subs_vals(dft_plus_u_section, &
1974  subsection_name="ENFORCE_OCCUPATION")
1975  subsection_enabled = .false.
1976  CALL section_vals_val_get(enforce_occupation_section, &
1977  keyword_name="_SECTION_PARAMETERS_", &
1978  l_val=subsection_enabled)
1979  IF (subsection_enabled) THEN
1980  NULLIFY (nelec)
1981  CALL section_vals_val_get(enforce_occupation_section, &
1982  keyword_name="NELEC", &
1983  r_vals=nelec)
1984  nspin = SIZE(nelec)
1985  ALLOCATE (qs_kind%dft_plus_u%nelec(nspin))
1986  qs_kind%dft_plus_u%nelec(:) = nelec(:)
1987  NULLIFY (orbitals)
1988  CALL section_vals_val_get(enforce_occupation_section, &
1989  keyword_name="ORBITALS", &
1990  i_vals=orbitals)
1991  norbitals = SIZE(orbitals)
1992  IF (norbitals <= 0 .OR. norbitals > 2*l + 1) &
1993  CALL cp_abort(__location__, "DFT+U| Invalid number of ORBITALS specified: "// &
1994  "1 to 2*L+1 integer numbers are expected")
1995  ALLOCATE (qs_kind%dft_plus_u%orbitals(norbitals))
1996  qs_kind%dft_plus_u%orbitals(:) = orbitals(:)
1997  NULLIFY (orbitals)
1998  DO m = 1, norbitals
1999  IF (qs_kind%dft_plus_u%orbitals(m) > l) &
2000  cpabort("DFT+U| Invalid orbital magnetic quantum number specified: m > l")
2001  IF (qs_kind%dft_plus_u%orbitals(m) < -l) &
2002  cpabort("DFT+U| Invalid orbital magnetic quantum number specified: m < -l")
2003  DO j = 1, norbitals
2004  IF (j /= m) THEN
2005  IF (qs_kind%dft_plus_u%orbitals(j) == qs_kind%dft_plus_u%orbitals(m)) &
2006  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  r_val=qs_kind%dft_plus_u%eps_scf)
2013  CALL section_vals_val_get(enforce_occupation_section, &
2014  keyword_name="MAX_SCF", &
2015  i_val=i)
2016  qs_kind%dft_plus_u%max_scf = max(-1, i)
2017  CALL section_vals_val_get(enforce_occupation_section, &
2018  keyword_name="SMEAR", &
2019  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  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  explicit_basis = .false.
2031  IF (k_rep > 0) THEN
2032  basis_section => section_vals_get_subs_vals(kind_section, "BASIS", i_rep_section=k_rep, &
2033  can_return_null=.true.)
2034  CALL section_vals_get(basis_section, explicit=explicit_basis)
2035  END IF
2036 
2037  explicit_potential = .false.
2038  IF (k_rep > 0) THEN
2039  potential_section => section_vals_get_subs_vals(kind_section, "POTENTIAL", &
2040  i_rep_section=k_rep, can_return_null=.true.)
2041  CALL section_vals_get(potential_section, explicit=explicit_potential)
2042  END IF
2043 
2044  explicit_kgpot = .false.
2045  IF (k_rep > 0) THEN
2046  kgpot_section => section_vals_get_subs_vals(kind_section, "KG_POTENTIAL", &
2047  i_rep_section=k_rep, can_return_null=.true.)
2048  CALL section_vals_get(kgpot_section, explicit=explicit_kgpot)
2049  END IF
2050 
2051  SELECT CASE (method_id)
2054  ! Allocate all_potential
2055  CALL allocate_potential(qs_kind%all_potential)
2056  CALL set_default_all_potential(qs_kind%all_potential, z, zeff_correction)
2057  CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
2058  IF (.NOT. ASSOCIATED(elec_conf)) THEN
2059  CALL get_potential(potential=qs_kind%all_potential, elec_conf=elec_conf)
2060  CALL set_qs_kind(qs_kind, elec_conf=elec_conf)
2061  END IF
2062  cpassert(.NOT. qs_kind%floating)
2063  IF (qs_kind%ghost) THEN
2064  CALL get_qs_kind(qs_kind=qs_kind, elec_conf=elec_conf)
2065  elec_conf(:) = 0
2066  CALL get_potential(potential=qs_kind%all_potential, &
2067  elec_conf=elec_conf)
2068  elec_conf(:) = 0
2069  CALL set_potential(potential=qs_kind%all_potential, &
2070  zeff=0.0_dp, &
2071  zeff_correction=0.0_dp)
2072  END IF
2073 
2074  ! Basis set (Parameters)
2075  ! Setup proper semiempirical parameters
2076  check = .NOT. ASSOCIATED(qs_kind%se_parameter)
2077  cpassert(check)
2078  CALL semi_empirical_create(qs_kind%se_parameter)
2079  ! Check if we allow p-orbitals on H
2080  SELECT CASE (z)
2081  CASE (1)
2082  IF (k_rep > 0) THEN
2083  CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
2084  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  CALL section_vals_val_get(dft_section, "QS%SE%STO_NG", i_val=ngauss)
2091  CALL se_param_set_default(qs_kind%se_parameter, z, method_id)
2092  NULLIFY (tmp_basis_set)
2093  CALL init_se_param(qs_kind%se_parameter, tmp_basis_set, ngauss)
2094  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  zeff=qs_kind%se_parameter%zeff, zeff_correction=zeff_correction)
2097  qs_kind%se_parameter%zeff = qs_kind%se_parameter%zeff - zeff_correction
2098 
2099  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  trim(qs_kind%name)//"> will be ignored!")
2104 
2105  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  trim(qs_kind%name)//"> will be ignored!")
2110 
2111  CASE (do_method_dftb)
2112  ! Allocate all_potential
2113  CALL allocate_potential(qs_kind%all_potential)
2114  CALL set_default_all_potential(qs_kind%all_potential, z, zeff_correction)
2115  CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
2116  IF (.NOT. ASSOCIATED(elec_conf)) THEN
2117  CALL get_potential(potential=qs_kind%all_potential, elec_conf=elec_conf)
2118  CALL set_qs_kind(qs_kind, elec_conf=elec_conf)
2119  END IF
2120  cpassert(.NOT. qs_kind%floating)
2121  IF (qs_kind%ghost) THEN
2122  CALL get_qs_kind(qs_kind=qs_kind, elec_conf=elec_conf)
2123  elec_conf(:) = 0
2124  CALL get_potential(potential=qs_kind%all_potential, &
2125  elec_conf=elec_conf)
2126  elec_conf(:) = 0
2127  CALL set_potential(potential=qs_kind%all_potential, &
2128  zeff=0.0_dp, &
2129  zeff_correction=0.0_dp)
2130  END IF
2131 
2132  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  trim(qs_kind%name)//"> will be ignored!")
2137 
2138  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  trim(qs_kind%name)//"> will be ignored!")
2143 
2144  CASE (do_method_xtb)
2145  ! Allocate all_potential
2146  CALL allocate_potential(qs_kind%all_potential)
2147  CALL set_default_all_potential(qs_kind%all_potential, z, zeff_correction)
2148  CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
2149  IF (.NOT. ASSOCIATED(elec_conf)) THEN
2150  CALL get_potential(potential=qs_kind%all_potential, elec_conf=elec_conf)
2151  CALL set_qs_kind(qs_kind, elec_conf=elec_conf)
2152  END IF
2153  cpassert(.NOT. qs_kind%floating)
2154  IF (qs_kind%ghost) THEN
2155  CALL get_qs_kind(qs_kind=qs_kind, elec_conf=elec_conf)
2156  elec_conf(:) = 0
2157  CALL get_potential(potential=qs_kind%all_potential, &
2158  elec_conf=elec_conf)
2159  elec_conf(:) = 0
2160  CALL set_potential(potential=qs_kind%all_potential, &
2161  zeff=0.0_dp, &
2162  zeff_correction=0.0_dp)
2163  END IF
2164 
2165  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  trim(qs_kind%name)//"> will be ignored!")
2170 
2171  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  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  IF (potential_name /= '') THEN
2181  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  "for the atomic kind <"//trim(qs_kind%name))
2187  CASE ("GTH")
2188  IF (potential_fn_kind == "-") THEN
2189  CALL section_vals_val_get(dft_section, "POTENTIAL_FILE_NAME", c_val=potential_file_name)
2190  ELSE
2191  potential_file_name = potential_fn_kind
2192  END IF
2193  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  potential_file_name, potential_section, update_input)
2197  CALL set_potential(qs_kind%gth_potential, z=z)
2198  CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
2199  IF (.NOT. ASSOCIATED(elec_conf)) THEN
2200  CALL get_potential(potential=qs_kind%gth_potential, elec_conf=elec_conf)
2201  CALL set_qs_kind(qs_kind, elec_conf=elec_conf)
2202  ELSE
2203  CALL set_potential(potential=qs_kind%gth_potential, elec_conf=elec_conf)
2204  END IF
2205  CASE ("UPF")
2206  ALLOCATE (qs_kind%upf_potential)
2207  CALL atom_read_upf(qs_kind%upf_potential, potential_name, read_header=.true.)
2208  CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
2209  IF (.NOT. ASSOCIATED(elec_conf)) THEN
2210  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  trim(qs_kind%name))
2218  END SELECT
2219  ELSE
2220  CALL cp_abort(__location__, &
2221  "No potential type was defined for the "// &
2222  "atomic kind <"//trim(qs_kind%name)//">")
2223  END IF
2224 
2225  CASE DEFAULT
2226 
2227  ! set ngauss for STO expansion
2228  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  DO i = 1, nb_rep
2232  SELECT CASE (basis_set_form(i))
2233  CASE ("GTO")
2234  NULLIFY (tmp_basis_set)
2235  CALL allocate_gto_basis_set(tmp_basis_set)
2236  CALL read_gto_basis_set(qs_kind%element_symbol, basis_set_name(i), &
2237  tmp_basis_set, para_env, dft_section)
2238  CASE ("STO")
2239  NULLIFY (sto_basis_set)
2240  CALL allocate_sto_basis_set(sto_basis_set)
2241  CALL read_sto_basis_set(qs_kind%element_symbol, basis_set_name(i), &
2242  sto_basis_set, para_env, dft_section)
2243  NULLIFY (tmp_basis_set)
2244  CALL create_gto_from_sto_basis(sto_basis_set, tmp_basis_set, ngauss)
2245  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  "for atomic kind <"//trim(qs_kind%name)//">")
2250  END SELECT
2251  tmp = basis_set_type(i)
2252  CALL uppercase(tmp)
2253  CALL add_basis_set_to_container(qs_kind%basis_sets, tmp_basis_set, tmp)
2254  END DO
2255  ! now explicit basis sets
2256  IF (explicit_basis) THEN
2257  CALL section_vals_get(basis_section, n_repetition=nexp)
2258  DO i = 1, nexp
2259  NULLIFY (tmp_basis_set)
2260  CALL allocate_gto_basis_set(tmp_basis_set)
2261  CALL read_gto_basis_set(qs_kind%element_symbol, basis_type, &
2262  tmp_basis_set, basis_section, i, dft_section)
2263  tmp = basis_type
2264  CALL uppercase(tmp)
2265  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  DO i = 1, SIZE(qs_kind%basis_sets)
2270  NULLIFY (tmp_basis_set)
2271  CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis_set, &
2272  inumbas=i, basis_type=basis_type)
2273  IF (basis_type == "") cycle
2274  jj = i
2275  DO j = i + 1, SIZE(qs_kind%basis_sets)
2276  jj = jj + 1
2277  NULLIFY (sup_basis_set)
2278  CALL get_basis_from_container(qs_kind%basis_sets, basis_set=sup_basis_set, &
2279  inumbas=jj, basis_type=tmp)
2280  IF (basis_type == tmp) THEN
2281  ! we found a match, combine the basis sets and delete the second
2282  CALL combine_basis_sets(tmp_basis_set, sup_basis_set)
2283  CALL remove_basis_from_container(qs_kind%basis_sets, jj)
2284  jj = jj - 1
2285  END IF
2286  END DO
2287  NULLIFY (sup_basis_set)
2288  END DO
2289 
2290  ! check that we have an orbital basis set
2291  nobasis = .true.
2292  DO i = 1, SIZE(qs_kind%basis_sets)
2293  NULLIFY (tmp_basis_set)
2294  CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis_set, &
2295  inumbas=i, basis_type=basis_type)
2296  IF (basis_type == "ORB") nobasis = .false.
2297  END DO
2298  IF (nobasis) THEN
2299  CALL cp_abort(__location__, &
2300  "No basis set type was defined for the "// &
2301  "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  IF (qs_kind%ghost .OR. qs_kind%floating) THEN
2306  IF (ASSOCIATED(qs_kind%elec_conf)) qs_kind%elec_conf = 0
2307  ELSE
2308  ! Allocate and initialise the potential data set structure
2309  IF ((potential_name /= '') .OR. explicit_potential) THEN
2310  ! determine the pseudopotential file to search
2311  IF (potential_fn_kind == "-") THEN
2312  CALL section_vals_val_get(dft_section, "POTENTIAL_FILE_NAME", c_val=potential_file_name)
2313  ELSE
2314  potential_file_name = potential_fn_kind
2315  END IF
2316  !
2317  SELECT CASE (trim(potential_type))
2318  CASE ("ALL")
2319  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  potential_file_name, potential_section, update_input)
2323  CALL set_potential(qs_kind%all_potential, z=z)
2324  CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
2325  IF (.NOT. ASSOCIATED(elec_conf)) THEN
2326  CALL get_potential(potential=qs_kind%all_potential, elec_conf=elec_conf)
2327  CALL set_qs_kind(qs_kind, elec_conf=elec_conf)
2328  ELSE
2329  CALL set_potential(potential=qs_kind%all_potential, elec_conf=elec_conf)
2330  END IF
2331  CASE ("GTH")
2332  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  potential_file_name, potential_section, update_input)
2336  CALL set_potential(qs_kind%gth_potential, z=z)
2337  CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
2338  IF (.NOT. ASSOCIATED(elec_conf)) THEN
2339  CALL get_potential(potential=qs_kind%gth_potential, elec_conf=elec_conf)
2340  CALL set_qs_kind(qs_kind, elec_conf=elec_conf)
2341  ELSE
2342  CALL set_potential(potential=qs_kind%gth_potential, elec_conf=elec_conf)
2343  END IF
2344  CASE ("ECP")
2345  CALL allocate_potential(qs_kind%sgp_potential)
2346  CALL get_potential(qs_kind%sgp_potential, description=description)
2347  CALL read_ecp_potential(ptable(z)%symbol, ecppot, &
2348  potential_name, potential_file_name, potential_section)
2349  IF (ecp_semi_local) THEN
2350  description(1) = "Semi-local Gaussian pseudopotential "
2351  description(2) = "ECP "//trim(potential_name)
2352  description(3) = "LIBGRPP: A. V. Oleynichenko et al., Symmetry 15 197 2023"
2353  description(4) = " "
2354  ELSE
2355  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  has_nlcc=.false.)
2361  CALL set_potential(qs_kind%sgp_potential, sl_lmax=ecppot%lmax, &
2362  npot=ecppot%npot, nrpot=ecppot%nrpot, apot=ecppot%apot, bpot=ecppot%bpot)
2363  ! convert PP
2364  IF (.NOT. ecp_semi_local) THEN
2365  cpabort("ECPs are only well tested in their semi-local form")
2366  CALL get_qs_kind(qs_kind, basis_set=orb_basis_set)
2367  CALL sgp_construction(sgp_pot=sgppot, ecp_pot=ecppot, orb_basis=orb_basis_set, error=error)
2368  IF (iounit > 0) THEN
2369  WRITE (iounit, "(/,T2,'PP Transformation for ',A)") trim(ecppot%pname)
2370  IF (sgppot%has_local) THEN
2371  WRITE (iounit, "(T8,'Accuracy for local part:',T41,F10.3,'%',T61,F20.12)") error(4), error(1)
2372  END IF
2373  IF (sgppot%has_nonlocal) THEN
2374  WRITE (iounit, "(T8,'Accuracy for nonlocal part:',T41,F10.3,'%',T61,F20.12)") error(5), error(2)
2375  END IF
2376  IF (sgppot%has_nlcc) THEN
2377  WRITE (iounit, "(T8,'Accuracy for NLCC density:',T61,F20.12)") error(3)
2378  END IF
2379  END IF
2380  END IF
2381  IF (sgppot%has_nonlocal) THEN
2382  CALL set_potential(qs_kind%sgp_potential, n_nonlocal=sgppot%n_nonlocal, lmax=sgppot%lmax, &
2383  is_nonlocal=sgppot%is_nonlocal)
2384  nnl = sgppot%n_nonlocal
2385  nppnl = 0
2386  DO l = 0, sgppot%lmax
2387  nppnl = nppnl + nnl*nco(l)
2388  END DO
2389  l = sgppot%lmax
2390  ALLOCATE (a_nl(nnl), h_nl(nnl, 0:l), c_nl(nnl, nnl, 0:l))
2391  a_nl(:) = sgppot%a_nonlocal(:)
2392  h_nl(:, :) = sgppot%h_nonlocal(:, :)
2393  DO l = 0, sgppot%lmax
2394  c_nl(:, :, l) = sgppot%c_nonlocal(:, :, l)*sqrt(2._dp*l + 1.0_dp)
2395  END DO
2396  CALL set_potential(qs_kind%sgp_potential, nppnl=nppnl, a_nonlocal=a_nl, h_nonlocal=h_nl, c_nonlocal=c_nl)
2397  ELSE
2398  CALL set_potential(qs_kind%sgp_potential, n_nonlocal=0, lmax=-1, is_nonlocal=sgppot%is_nonlocal)
2399  CALL set_potential(qs_kind%sgp_potential, nppnl=0)
2400  END IF
2401  !
2402  cpassert(.NOT. sgppot%has_local)
2403  cpassert(.NOT. sgppot%has_nlcc)
2404  ! core
2405  rc = 0.5_dp*qs_kind%covalent_radius*angstrom
2406  rc = max(rc, 0.2_dp)
2407  rc = min(rc, 1.0_dp)
2408  alpha = 1.0_dp/(2.0_dp*rc**2)
2409  ccore = ecppot%zion*sqrt((alpha/pi)**3)
2410  CALL set_potential(qs_kind%sgp_potential, alpha_core_charge=alpha, ccore_charge=ccore, &
2411  core_charge_radius=rc)
2412  CALL atom_sgp_release(sgppot)
2413  CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
2414  IF (.NOT. ASSOCIATED(elec_conf)) THEN
2415  CALL set_qs_kind(qs_kind, elec_conf=ecppot%econf)
2416  END IF
2417  CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
2418  CALL set_potential(qs_kind%sgp_potential, elec_conf=elec_conf)
2419  CASE ("UPF")
2420  CALL allocate_potential(qs_kind%sgp_potential)
2421  CALL get_potential(qs_kind%sgp_potential, description=description)
2422  description(4) = "UPF "//trim(potential_name)
2423  CALL atom_read_upf(upfpot, potential_name)
2424  CALL set_potential(qs_kind%sgp_potential, name=upfpot%pname, description=description, &
2425  zeff=upfpot%zion, z=z, has_nlcc=upfpot%core_correction)
2426  ! convert pp
2427  CALL sgp_construction(sgp_pot=sgppot, upf_pot=upfpot, error=error)
2428  IF (iounit > 0) THEN
2429  WRITE (iounit, "(/,T2,'PP Transformation for ',A)") trim(upfpot%pname)
2430  IF (sgppot%has_local) THEN
2431  WRITE (iounit, "(T8,'Accuracy for local part:',T61,F20.12)") error(1)
2432  END IF
2433  IF (sgppot%has_nonlocal) THEN
2434  WRITE (iounit, "(T8,'Accuracy for nonlocal part:',T61,F20.12)") error(2)
2435  END IF
2436  IF (sgppot%has_nlcc) THEN
2437  WRITE (iounit, "(T8,'Accuracy for NLCC density:',T61,F20.12)") error(3)
2438  END IF
2439  END IF
2440  IF (sgppot%has_nonlocal) THEN
2441  CALL set_potential(qs_kind%sgp_potential, n_nonlocal=sgppot%n_nonlocal, lmax=sgppot%lmax, &
2442  is_nonlocal=sgppot%is_nonlocal)
2443  nnl = sgppot%n_nonlocal
2444  nppnl = 0
2445  DO l = 0, sgppot%lmax
2446  nppnl = nppnl + nnl*nco(l)
2447  END DO
2448  l = sgppot%lmax
2449  ALLOCATE (a_nl(nnl), h_nl(nnl, 0:l), c_nl(nnl, nnl, 0:l))
2450  a_nl(:) = sgppot%a_nonlocal(:)
2451  h_nl(:, :) = sgppot%h_nonlocal(:, :)
2452  c_nl(:, :, :) = sgppot%c_nonlocal(:, :, :)
2453  CALL set_potential(qs_kind%sgp_potential, nppnl=nppnl, a_nonlocal=a_nl, h_nonlocal=h_nl, c_nonlocal=c_nl)
2454  ELSE
2455  CALL set_potential(qs_kind%sgp_potential, n_nonlocal=0, lmax=-1, is_nonlocal=sgppot%is_nonlocal)
2456  CALL set_potential(qs_kind%sgp_potential, nppnl=0)
2457  END IF
2458  cpassert(sgppot%has_local)
2459  ! core
2460  rc = sgppot%ac_local
2461  alpha = 1.0_dp/(2.0_dp*rc**2)
2462  ccore = upfpot%zion*sqrt((alpha/pi)**3)
2463  CALL set_potential(qs_kind%sgp_potential, alpha_core_charge=alpha, ccore_charge=ccore, &
2464  core_charge_radius=rc)
2465  ! local potential
2466  nloc = sgppot%n_local
2467  ALLOCATE (aloc(nloc), cloc(nloc))
2468  aloc(1:nloc) = sgppot%a_local(1:nloc)
2469  cloc(1:nloc) = sgppot%c_local(1:nloc)
2470  CALL set_potential(qs_kind%sgp_potential, n_local=nloc, a_local=aloc, c_local=cloc)
2471  IF (sgppot%has_nlcc) THEN
2472  nlcc = sgppot%n_nlcc
2473  ALLOCATE (anlcc(nlcc), cnlcc(nlcc))
2474  anlcc(1:nlcc) = sgppot%a_nlcc(1:nlcc)
2475  cnlcc(1:nlcc) = sgppot%c_nlcc(1:nlcc)
2476  CALL set_potential(qs_kind%sgp_potential, has_nlcc=.true., n_nlcc=nlcc, a_nlcc=anlcc, c_nlcc=cnlcc)
2477  END IF
2478  CALL set_potential(qs_kind%sgp_potential, z=z)
2479  CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
2480  IF (.NOT. ASSOCIATED(elec_conf)) THEN
2481  CALL set_qs_kind(qs_kind, elec_conf=upfpot%econf)
2482  END IF
2483  CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
2484  CALL set_potential(qs_kind%sgp_potential, elec_conf=elec_conf)
2485  CALL atom_release_upf(upfpot)
2486  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  trim(qs_kind%name))
2493  END SELECT
2494  ELSE
2495  CALL cp_abort(__location__, &
2496  "No potential type was defined for the "// &
2497  "atomic kind <"//trim(qs_kind%name)//">")
2498  END IF
2499 
2500  CALL check_potential_basis_compatibility(qs_kind)
2501 
2502  ! Allocate and initialise the potential data set structure
2503  IF ((kgpot_name /= '') .OR. explicit_kgpot) THEN
2504  ipos = index(kgpot_name, "-")
2505  IF (ipos > 1) THEN
2506  kgpot_type = kgpot_name(:ipos - 1)
2507  ELSE
2508  kgpot_type = kgpot_name
2509  END IF
2510  CALL uppercase(kgpot_type)
2511 
2512  SELECT CASE (trim(kgpot_type))
2513  CASE ("TNADD")
2514  ! determine the pseudopotential file to search
2515  IF (kg_potential_fn_kind == "-") THEN
2516  CALL section_vals_val_get(dft_section, "POTENTIAL_FILE_NAME", c_val=potential_file_name)
2517  ELSE
2518  potential_file_name = kg_potential_fn_kind
2519  END IF
2520  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  potential_file_name, kgpot_section, update_input)
2524  CASE ("NONE")
2525  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  trim(qs_kind%name))
2532  END SELECT
2533  END IF
2534  END IF
2535  END SELECT
2536 
2537  CALL timestop(handle)
2538 
2539  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  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  CALL get_qs_kind(qs_kind, name=name, gth_potential=gth_potential, basis_set=basis_set)
2555 
2556  npp = -1; nbs = -1
2557  IF (ASSOCIATED(gth_potential)) &
2558  npp = parse_valence_electrons(gth_potential%aliases)
2559  IF (ASSOCIATED(basis_set)) &
2560  nbs = parse_valence_electrons(basis_set%aliases)
2561 
2562  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  " were optimized for different valence electron numbers.")
2565 
2566  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  FUNCTION parse_valence_electrons(string) RESULT(n)
2575  CHARACTER(*) :: string
2576  INTEGER :: n
2577 
2578  INTEGER :: i, istat, j
2579 
2580  i = index(string, "-Q", .true.)
2581  IF (i == 0) THEN
2582  n = -1
2583  ELSE
2584  j = scan(string(i + 2:), "- ")
2585  READ (string(i + 2:i + j), '(I3)', iostat=istat) n
2586  IF (istat /= 0) n = -1
2587  END IF
2588 
2589  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  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  CALL timeset(routinen, handle)
2613 
2614  IF (ASSOCIATED(qs_kind_set)) cpabort("create_qs_kind_set: qs_kind_set already associated")
2615  IF (.NOT. ASSOCIATED(atomic_kind_set)) cpabort("create_qs_kind_set: atomic_kind_set not associated")
2616 
2617  no_fail = .false.
2618 
2619  ! Between all methods only SE and DFTB/xTB may not need a KIND section.
2620  CALL section_vals_val_get(force_env_section, "METHOD", i_val=method)
2621  IF (method == do_qs) THEN
2622  CALL section_vals_val_get(force_env_section, "DFT%QS%METHOD", i_val=qs_method)
2623  SELECT CASE (qs_method)
2626  no_fail = .true.
2627  CASE (do_method_dftb)
2628  no_fail = .true.
2629  CASE (do_method_xtb)
2630  no_fail = .true.
2631  END SELECT
2632  ELSE IF (method == do_sirius) THEN
2633  qs_method = do_method_pw
2634  ELSE
2635  qs_method = method
2636  END IF
2637 
2638  nkind = SIZE(atomic_kind_set)
2639  ALLOCATE (qs_kind_set(nkind))
2640 
2641  DO ikind = 1, nkind
2642  qs_kind_set(ikind)%name = atomic_kind_set(ikind)%name
2643  qs_kind_set(ikind)%element_symbol = atomic_kind_set(ikind)%element_symbol
2644  qs_kind_set(ikind)%natom = atomic_kind_set(ikind)%natom
2645  CALL read_qs_kind(qs_kind_set(ikind), kind_section, para_env, force_env_section, no_fail, qs_method)
2646  END DO
2647 
2648  CALL timestop(handle)
2649 
2650  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  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  IF (dft_control%qs_control%semi_empirical) THEN
2671  CALL get_qs_kind(qs_kind, se_parameter=se_parameter)
2672  cpassert(ASSOCIATED(se_parameter))
2673  CALL get_se_param(se_parameter, defined=defined)
2674  cpassert(defined)
2675  CALL write_se_param(se_parameter, subsys_section)
2676  ELSE IF (dft_control%qs_control%dftb) THEN
2677  CALL get_qs_kind(qs_kind, dftb_parameter=dftb_parameter)
2678  cpassert(ASSOCIATED(dftb_parameter))
2679  CALL get_dftb_atom_param(dftb_parameter, defined=defined)
2680  cpassert(defined)
2681  CALL write_dftb_atom_param(dftb_parameter, subsys_section)
2682  ELSE IF (dft_control%qs_control%xtb) THEN
2683  CALL get_qs_kind(qs_kind, xtb_parameter=xtb_parameter)
2684  cpassert(ASSOCIATED(xtb_parameter))
2685  CALL write_xtb_atom_param(xtb_parameter, subsys_section)
2686  END IF
2687 
2688  END SUBROUTINE check_qs_kind
2689 
2690 ! **************************************************************************************************
2691 !> \brief ...
2692 !> \param qs_kind_set ...
2693 !> \param dft_control ...
2694 !> \param subsys_section ...
2695 ! **************************************************************************************************
2696  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  CALL timeset(routinen, handle)
2708  IF (ASSOCIATED(qs_kind_set)) THEN
2709  nkind = SIZE(qs_kind_set)
2710  DO ikind = 1, nkind
2711  qs_kind => qs_kind_set(ikind)
2712  CALL check_qs_kind(qs_kind, dft_control, subsys_section)
2713  END DO
2714  IF (dft_control%qs_control%xtb) THEN
2715  CALL write_xtb_kab_param(qs_kind_set, subsys_section, &
2716  dft_control%qs_control%xtb_control)
2717  END IF
2718  ELSE
2719  cpabort("The pointer qs_kind_set is not associated")
2720  END IF
2721  CALL timestop(handle)
2722  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  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  NULLIFY (logger)
2743  logger => cp_get_default_logger()
2744  IF (btest(cp_print_key_should_output(logger%iter_info, subsys_section, &
2745  "PRINT%KINDS/POTENTIAL"), cp_p_file)) THEN
2746 
2747  io_unit = cp_print_key_unit_nr(logger, subsys_section, "PRINT%KINDS", extension=".Log")
2748  IF (io_unit > 0) THEN
2749 
2750  WRITE (io_unit, "(/,T2,A)") "xTB| Kab parameters"
2751  nkind = SIZE(qs_kind_set)
2752  DO ikind = 1, nkind
2753  qs_kinda => qs_kind_set(ikind)
2754  CALL get_qs_kind(qs_kinda, xtb_parameter=xtb_parameter_a)
2755  CALL get_xtb_atom_param(xtb_parameter_a, aname=aname, z=za)
2756  DO jkind = ikind, nkind
2757  qs_kindb => qs_kind_set(jkind)
2758  CALL get_qs_kind(qs_kindb, xtb_parameter=xtb_parameter_b)
2759  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  " Kab:", trim(aname), trim(bname), xtb_set_kab(za, zb, xtb_control)
2762  END DO
2763  END DO
2764  WRITE (io_unit, *)
2765 
2766  END IF
2767 
2768  CALL cp_print_key_finished_output(io_unit, logger, subsys_section, "PRINT%KINDS")
2769  END IF
2770 
2771  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  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  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  IF (PRESENT(dftb_parameter)) qs_kind%dftb_parameter => dftb_parameter
2816  IF (PRESENT(xtb_parameter)) qs_kind%xtb_parameter => xtb_parameter
2817  IF (PRESENT(elec_conf)) THEN
2818  IF (ASSOCIATED(qs_kind%elec_conf)) THEN
2819  DEALLOCATE (qs_kind%elec_conf)
2820  END IF
2821  ALLOCATE (qs_kind%elec_conf(0:SIZE(elec_conf) - 1))
2822  qs_kind%elec_conf(:) = elec_conf(:)
2823  END IF
2824  IF (PRESENT(paw_atom)) qs_kind%paw_atom = paw_atom
2825  IF (PRESENT(hard_radius)) qs_kind%hard_radius = hard_radius
2826  IF (PRESENT(hard0_radius)) qs_kind%hard0_radius = hard0_radius
2827  IF (PRESENT(covalent_radius)) qs_kind%covalent_radius = covalent_radius
2828  IF (PRESENT(vdw_radius)) qs_kind%vdw_radius = vdw_radius
2829  IF (PRESENT(lmax_rho0)) qs_kind%lmax_rho0 = lmax_rho0
2830  IF (PRESENT(zeff)) THEN
2831  IF (ASSOCIATED(qs_kind%all_potential)) THEN
2832  CALL set_potential(potential=qs_kind%all_potential, zeff=zeff)
2833  ELSE IF (ASSOCIATED(qs_kind%gth_potential)) THEN
2834  CALL set_potential(potential=qs_kind%gth_potential, zeff=zeff)
2835  ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
2836  CALL set_potential(potential=qs_kind%sgp_potential, zeff=zeff)
2837  END IF
2838  END IF
2839  IF (PRESENT(ghost)) qs_kind%ghost = ghost
2840 
2841  IF (PRESENT(floating)) qs_kind%floating = floating
2842 
2843  IF (PRESENT(no_optimize)) qs_kind%no_optimize = no_optimize
2844 
2845  IF (PRESENT(dispersion)) qs_kind%dispersion => dispersion
2846 
2847  IF (PRESENT(u_minus_j)) THEN
2848  IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
2849  qs_kind%dft_plus_u%u_minus_j = u_minus_j
2850  END IF
2851  END IF
2852 
2853  IF (PRESENT(reltmat)) qs_kind%reltmat => reltmat
2854 
2855  IF (PRESENT(pao_basis_size)) qs_kind%pao_basis_size = pao_basis_size
2856 
2857  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  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  IF (output_unit > 0) THEN
2879 
2880  IF (ASSOCIATED(qs_kind)) THEN
2881  WRITE (unit=output_unit, fmt="(/,T2,I2,A,T57,A,T75,I6)") &
2882  kind_number, ". Atomic kind: "//trim(qs_kind%name), &
2883  "Number of atoms: ", qs_kind%natom
2884 
2885  DO ibas = 1, SIZE(qs_kind%basis_sets, 1)
2886  NULLIFY (tmp_basis)
2887  CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis, &
2888  inumbas=ibas, basis_type=basis_type)
2889  do_print = .true.
2890  SELECT CASE (basis_type)
2891  CASE DEFAULT
2892  bstring = "Basis Set"
2893  do_print = .false.
2894  CASE ("ORB")
2895  bstring = "Orbital Basis Set"
2896  CASE ("ORB_SOFT")
2897  bstring = "GAPW Soft Basis Set"
2898  do_print = .false.
2899  CASE ("AUX")
2900  bstring = "Auxiliary Basis Set"
2901  CASE ("MIN")
2902  bstring = "Minimal Basis Set"
2903  CASE ("RI_AUX")
2904  bstring = "RI Auxiliary Basis Set"
2905  CASE ("AUX_FIT")
2906  bstring = "Auxiliary Fit Basis Set"
2907  CASE ("LRI_AUX")
2908  bstring = "LRI Basis Set"
2909  CASE ("P_LRI_AUX")
2910  bstring = "LRI Basis Set for TDDFPT"
2911  CASE ("RI_XAS")
2912  bstring = "RI XAS Basis Set"
2913  CASE ("RI_HFX")
2914  bstring = "RI HFX Basis Set"
2915  END SELECT
2916 
2917  IF (do_print) THEN
2918  CALL write_orb_basis_set(tmp_basis, output_unit, bstring)
2919  END IF
2920 
2921  END DO
2922 
2923  IF (qs_kind%ghost) THEN
2924  WRITE (unit=output_unit, fmt="(/,T6,A)") &
2925  "The atoms of this atomic kind are GHOST atoms!"
2926  END IF
2927  IF (qs_kind%floating) THEN
2928  WRITE (unit=output_unit, fmt="(/,T6,A)") &
2929  "The atoms of this atomic kind are FLOATING BASIS FUNCTIONS."
2930  END IF
2931  IF (qs_kind%covalent_radius > 0.0_dp) THEN
2932  WRITE (unit=output_unit, fmt="(/,T8,A,T71,F10.3)") &
2933  "Atomic covalent radius [Angstrom]:", &
2934  qs_kind%covalent_radius*angstrom
2935  END IF
2936  IF (qs_kind%vdw_radius > 0.0_dp) THEN
2937  WRITE (unit=output_unit, fmt="(/,T8,A,T71,F10.3)") &
2938  "Atomic van der Waals radius [Angstrom]:", &
2939  qs_kind%vdw_radius*angstrom
2940  END IF
2941  IF (qs_kind%paw_atom) THEN
2942  WRITE (unit=output_unit, fmt="(/,T6,A)") &
2943  "The atoms of this atomic kind are PAW atoms (GAPW):"
2944  WRITE (unit=output_unit, fmt="(T8,A,T71,F10.3)") &
2945  "Hard Gaussian function radius:", qs_kind%hard_radius, &
2946  "Rho0 radius:", qs_kind%hard0_radius, &
2947  "Maximum GTO radius used for PAW projector construction:", &
2948  qs_kind%max_rad_local
2949  NULLIFY (tmp_basis)
2950  CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis, &
2951  basis_type="ORB_SOFT")
2952  CALL write_orb_basis_set(tmp_basis, output_unit, "GAPW Soft Basis Set")
2953  END IF
2954  ! Potentials
2955  IF (ASSOCIATED(qs_kind%all_potential)) CALL write_potential(qs_kind%all_potential, output_unit)
2956  IF (ASSOCIATED(qs_kind%gth_potential)) CALL write_potential(qs_kind%gth_potential, output_unit)
2957  IF (ASSOCIATED(qs_kind%sgp_potential)) CALL write_potential(qs_kind%sgp_potential, output_unit)
2958  IF (ASSOCIATED(qs_kind%tnadd_potential)) CALL write_potential(qs_kind%tnadd_potential, output_unit)
2959  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  "A DFT+U correction is applied to atoms of this atomic kind:", &
2962  "Angular quantum momentum number L:", qs_kind%dft_plus_u%l, &
2963  "U(eff) = (U - J) value in [eV]:", qs_kind%dft_plus_u%u_minus_j_target*evolt
2964  IF (qs_kind%dft_plus_u%u_ramping > 0.0_dp) THEN
2965  IF (qs_kind%dft_plus_u%init_u_ramping_each_scf) THEN
2966  yon = "YES"
2967  ELSE
2968  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  "Increment for U ramping in [eV]:", qs_kind%dft_plus_u%u_ramping*evolt, &
2972  "SCF threshold value for U ramping:", qs_kind%dft_plus_u%eps_u_ramping, &
2973  "Set U ramping value to zero before each wavefunction optimisation:", yon
2974  END IF
2975  IF (ASSOCIATED(qs_kind%dft_plus_u%orbitals)) THEN
2976  WRITE (unit=output_unit, fmt="(T8,A)") &
2977  "An initial orbital occupation is requested:"
2978  IF (ASSOCIATED(qs_kind%dft_plus_u%nelec)) THEN
2979  IF (any(qs_kind%dft_plus_u%nelec(:) >= 0.5_dp)) THEN
2980  IF (SIZE(qs_kind%dft_plus_u%nelec) > 1) THEN
2981  WRITE (unit=output_unit, fmt="(T9,A,T75,F6.2)") &
2982  "Number of alpha electrons:", &
2983  qs_kind%dft_plus_u%nelec(1), &
2984  "Number of beta electrons:", &
2985  qs_kind%dft_plus_u%nelec(2)
2986  ELSE
2987  WRITE (unit=output_unit, fmt="(T9,A,T75,F6.2)") &
2988  "Number of electrons:", &
2989  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  "Preferred (initial) orbital occupation order (orbital M values):", &
2995  qs_kind%dft_plus_u%orbitals(:)
2996  WRITE (unit=output_unit, fmt="(T9,A,T71,ES10.3,/,T9,A,T76,I5)") &
2997  "Threshold value for the SCF convergence criterion:", &
2998  qs_kind%dft_plus_u%eps_scf, &
2999  "Number of initial SCF iterations:", &
3000  qs_kind%dft_plus_u%max_scf
3001  IF (qs_kind%dft_plus_u%smear) THEN
3002  WRITE (unit=output_unit, fmt="(T9,A)") &
3003  "A smearing of the orbital occupations will be performed"
3004  END IF
3005  END IF
3006  END IF
3007  ELSE
3008  cpabort("")
3009  END IF
3010 
3011  END IF
3012 
3013  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  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  CALL timeset(routinen, handle)
3033 
3034  NULLIFY (logger)
3035  logger => cp_get_default_logger()
3036  output_unit = cp_print_key_unit_nr(logger, subsys_section, &
3037  "PRINT%KINDS", extension=".Log")
3038  IF (output_unit > 0) THEN
3039  IF (ASSOCIATED(qs_kind_set)) THEN
3040  WRITE (unit=output_unit, fmt="(/,/,T2,A)") "ATOMIC KIND INFORMATION"
3041  nkind = SIZE(qs_kind_set)
3042  DO ikind = 1, nkind
3043  qs_kind => qs_kind_set(ikind)
3044  CALL write_qs_kind(qs_kind, ikind, output_unit)
3045  END DO
3046  ELSE
3047  cpabort("")
3048  END IF
3049  END IF
3050 
3051  CALL cp_print_key_finished_output(output_unit, logger, subsys_section, &
3052  "PRINT%KINDS")
3053 
3054  CALL timestop(handle)
3055 
3056  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  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  CALL timeset(routinen, handle)
3081 
3082  NULLIFY (logger)
3083  logger => cp_get_default_logger()
3084  output_unit = cp_print_key_unit_nr(logger, subsys_section, &
3085  "PRINT%KINDS/BASIS_SET", &
3086  extension=".Log")
3087  IF (output_unit > 0) THEN
3088  IF (ASSOCIATED(qs_kind_set)) THEN
3089  WRITE (unit=output_unit, fmt="(/,/,T2,A)") &
3090  "BASIS SET INFORMATION (Unnormalised Gaussian-type functions)"
3091  nkind = SIZE(qs_kind_set)
3092  DO ikind = 1, nkind
3093  qs_kind => qs_kind_set(ikind)
3094  WRITE (unit=output_unit, fmt="(/,T2,I2,A)") &
3095  ikind, ". Atomic kind: "//trim(qs_kind%name)
3096 
3097  DO ibas = 1, SIZE(qs_kind%basis_sets, 1)
3098  NULLIFY (tmp_basis)
3099  CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis, &
3100  inumbas=ibas, basis_type=basis_type)
3101  IF (basis_type == "") cycle
3102  SELECT CASE (basis_type)
3103  CASE DEFAULT
3104  bstring = "Basis Set"
3105  CASE ("ORB")
3106  bstring = "Orbital Basis Set"
3107  CASE ("ORB_SOFT")
3108  bstring = "GAPW Soft Basis Set"
3109  CASE ("AUX")
3110  bstring = "Auxiliary Basis Set"
3111  CASE ("MIN")
3112  bstring = "Minimal Basis Set"
3113  CASE ("RI_AUX")
3114  bstring = "RI Auxiliary Basis Set"
3115  CASE ("AUX_FIT")
3116  bstring = "Auxiliary Fit Basis Set"
3117  CASE ("LRI_AUX")
3118  bstring = "LRI Basis Set"
3119  CASE ("P_LRI_AUX")
3120  bstring = "LRI Basis Set for TDDFPT"
3121  CASE ("RI_HFX")
3122  bstring = "RI HFX Basis Set"
3123  END SELECT
3124 
3125  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  cpabort("")
3132  END IF
3133  END IF
3134 
3135  CALL cp_print_key_finished_output(output_unit, logger, subsys_section, &
3136  "PRINT%KINDS/BASIS_SET")
3137 
3138  CALL timestop(handle)
3139 
3140  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  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  INTEGER, DIMENSION(:), POINTER :: econf
3161  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  CALL get_atomic_kind(atomic_kind, z=z)
3168  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  addel=addel, laddel=laddel, naddel=naddel)
3175 
3176  ! electronic state
3177  nelem = 0
3178  ncore = 0
3179  ncalc = 0
3180  edelta = 0.0_dp
3181  IF (ASSOCIATED(gth_potential)) THEN
3182  CALL get_potential(gth_potential, elec_conf=econf)
3183  CALL set_pseudo_state(econf, z, ncalc, ncore, nelem)
3184  ELSE IF (ASSOCIATED(sgp_potential)) THEN
3185  CALL get_potential(sgp_potential, elec_conf=econf)
3186  CALL set_pseudo_state(econf, z, ncalc, ncore, nelem)
3187  ELSE
3188  DO l = 0, min(lmat, ubound(ptable(z)%e_conv, 1))
3189  ll = 2*(2*l + 1)
3190  nn = ptable(z)%e_conv(l)
3191  ii = 0
3192  DO
3193  ii = ii + 1
3194  IF (nn <= ll) THEN
3195  nelem(l, ii) = nn
3196  EXIT
3197  ELSE
3198  nelem(l, ii) = ll
3199  nn = nn - ll
3200  END IF
3201  END DO
3202  END DO
3203  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  IF (bs_occupation) THEN
3210  DO is = 1, 2
3211  DO i = 1, SIZE(addel, 1)
3212  ne = addel(i, is)
3213  l = laddel(i, is)
3214  nn = naddel(i, is) - l
3215  IF (ne /= 0) THEN
3216  IF (nn == 0) THEN
3217  DO ii = SIZE(nelem, 2), 1, -1
3218  IF (ncalc(l, ii) > 0) THEN
3219  IF ((ncalc(l, ii) + ne) < 2*(2*l + 1) + 1) THEN
3220  edelta(l, ii, is) = edelta(l, ii, is) + ne
3221  nn = ii
3222  ELSE
3223  edelta(l, ii + 1, is) = edelta(l, ii + 1, is) + ne
3224  nn = ii + 1
3225  END IF
3226  EXIT
3227  ELSE IF (ii == 1) THEN
3228  edelta(l, ii, is) = edelta(l, ii, is) + ne
3229  nn = ii
3230  END IF
3231  END DO
3232  ELSE
3233  edelta(l, nn, is) = edelta(l, nn, is) + ne
3234  END IF
3235  IF (ncalc(l, nn) + edelta(l, nn, is) < 0) THEN
3236  edelta(l, nn, is) = -ncalc(l, nn)
3237  END IF
3238  END IF
3239  END DO
3240  END DO
3241  edelta = 0.5_dp*edelta
3242  ELSE IF (magnetization /= 0.0_dp) THEN
3243  dmag = 0.5_dp*abs(magnetization)
3244  DO l = 0, min(lmat, ubound(ptable(z)%e_conv, 1))
3245  ll = 2*(2*l + 1)
3246  ii = 0
3247  DO i = 1, SIZE(ncalc, 2)
3248  IF (ncalc(l, i) == 0) cycle
3249  IF (ncalc(l, i) == ll) cycle
3250  IF (ncalc(l, i) > dmag .AND. (ll - ncalc(l, i)) > dmag) THEN
3251  ii = i
3252  EXIT
3253  END IF
3254  END DO
3255  IF (ii /= 0) THEN
3256  edelta(l, ii, 1) = magnetization*0.5_dp
3257  edelta(l, ii, 2) = -magnetization*0.5_dp
3258  EXIT
3259  END IF
3260  END DO
3261  IF (ii == 0) THEN
3262  CALL cp_abort(__location__, &
3263  "Magnetization value cannot be imposed for this atom type")
3264  END IF
3265  END IF
3266 
3267  IF (qs_kind%ghost .OR. qs_kind%floating) THEN
3268  nelem = 0
3269  ncore = 0
3270  ncalc = 0
3271  edelta = 0.0_dp
3272  END IF
3273 
3274  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  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  NULLIFY (logger)
3295  logger => cp_get_default_logger()
3296  iounit = cp_logger_get_default_io_unit(logger)
3297 
3298  econfx = 0
3299  econfx(0:SIZE(econf) - 1) = econf
3300  IF (sum(econf) >= 0) THEN
3301  lmin = min(lmat, ubound(ptable(z)%e_conv, 1))
3302  ! number of core electrons
3303  nc = z - sum(econf)
3304  ! setup ncore
3305  ncore = 0
3306  SELECT CASE (nc)
3307  CASE (0)
3308  CASE (2)
3309  ncore(0, 1) = 2
3310  CASE (10)
3311  ncore(0, 1) = 2
3312  ncore(0, 2) = 2
3313  ncore(1, 1) = 6
3314  CASE (18)
3315  ncore(0, 1) = 2
3316  ncore(0, 2) = 2
3317  ncore(0, 3) = 2
3318  ncore(1, 1) = 6
3319  ncore(1, 2) = 6
3320  CASE (28)
3321  ncore(0, 1) = 2
3322  ncore(0, 2) = 2
3323  ncore(0, 3) = 2
3324  ncore(1, 1) = 6
3325  ncore(1, 2) = 6
3326  ncore(2, 1) = 10
3327  CASE (36)
3328  ncore(0, 1) = 2
3329  ncore(0, 2) = 2
3330  ncore(0, 3) = 2
3331  ncore(0, 4) = 2
3332  ncore(1, 1) = 6
3333  ncore(1, 2) = 6
3334  ncore(1, 3) = 6
3335  ncore(2, 1) = 10
3336  CASE (46)
3337  ncore(0, 1) = 2
3338  ncore(0, 2) = 2
3339  ncore(0, 3) = 2
3340  ncore(0, 4) = 2
3341  ncore(1, 1) = 6
3342  ncore(1, 2) = 6
3343  ncore(1, 3) = 6
3344  ncore(2, 1) = 10
3345  ncore(2, 2) = 10
3346  CASE (54)
3347  ncore(0, 1) = 2
3348  ncore(0, 2) = 2
3349  ncore(0, 3) = 2
3350  ncore(0, 4) = 2
3351  ncore(0, 5) = 2
3352  ncore(1, 1) = 6
3353  ncore(1, 2) = 6
3354  ncore(1, 3) = 6
3355  ncore(1, 4) = 6
3356  ncore(2, 1) = 10
3357  ncore(2, 2) = 10
3358  CASE (60)
3359  ncore(0, 1) = 2
3360  ncore(0, 2) = 2
3361  ncore(0, 3) = 2
3362  ncore(0, 4) = 2
3363  ncore(1, 1) = 6
3364  ncore(1, 2) = 6
3365  ncore(1, 3) = 6
3366  ncore(2, 1) = 10
3367  ncore(2, 2) = 10
3368  ncore(3, 1) = 14
3369  CASE (68)
3370  ncore(0, 1) = 2
3371  ncore(0, 2) = 2
3372  ncore(0, 3) = 2
3373  ncore(0, 4) = 2
3374  ncore(0, 5) = 2
3375  ncore(1, 1) = 6
3376  ncore(1, 2) = 6
3377  ncore(1, 3) = 6
3378  ncore(1, 4) = 6
3379  ncore(2, 1) = 10
3380  ncore(2, 2) = 10
3381  ncore(3, 1) = 14
3382  CASE (78)
3383  ncore(0, 1) = 2
3384  ncore(0, 2) = 2
3385  ncore(0, 3) = 2
3386  ncore(0, 4) = 2
3387  ncore(0, 5) = 2
3388  ncore(1, 1) = 6
3389  ncore(1, 2) = 6
3390  ncore(1, 3) = 6
3391  ncore(1, 4) = 6
3392  ncore(2, 1) = 10
3393  ncore(2, 2) = 10
3394  ncore(2, 3) = 10
3395  ncore(3, 1) = 14
3396  CASE DEFAULT
3397  ncore(0, 1) = -1
3398  END SELECT
3399  ! special cases of double assignments
3400  IF (z == 65 .AND. econfx(3) == 0) THEN
3401  ! 4f in core for Tb
3402  ncore = 0
3403  ncore(0, 1) = -1
3404  END IF
3405  ! if there is still no core, check for special cases
3406  IF (ncore(0, 1) <= 0) THEN
3407  IF (z >= 58 .AND. z <= 71) THEN
3408  ! 4f-in-core PPs for lanthanides
3409  nc = z - sum(econf)
3410  ! setup ncore
3411  ncore = 0
3412  SELECT CASE (nc)
3413  CASE (29:42)
3414  ncore(0, 1) = 2
3415  ncore(0, 2) = 2
3416  ncore(0, 3) = 2
3417  ncore(1, 1) = 6
3418  ncore(1, 2) = 6
3419  ncore(2, 1) = 10
3420  ncore(3, 1) = nc - 28
3421  message = "A small-core pseudopotential with 4f-in-core is used for the lanthanide "// &
3422  trim(ptable(z)%symbol)
3423  cphint(trim(message))
3424  CASE (47:60)
3425  ncore(0, 1) = 2
3426  ncore(0, 2) = 2
3427  ncore(0, 3) = 2
3428  ncore(0, 4) = 2
3429  ncore(1, 1) = 6
3430  ncore(1, 2) = 6
3431  ncore(1, 3) = 6
3432  ncore(2, 1) = 10
3433  ncore(2, 2) = 10
3434  ncore(3, 1) = nc - 46
3435  message = "A medium-core pseudopotential with 4f-in-core is used for the lanthanide "// &
3436  trim(ptable(z)%symbol)
3437  cphint(trim(message))
3438  CASE DEFAULT
3439  ncore(0, 1) = -1
3440  END SELECT
3441  END IF
3442  END IF
3443  ! if the core is established, finish the setup
3444  IF (ncore(0, 1) >= 0) THEN
3445  DO l = 0, lmin
3446  ll = 2*(2*l + 1)
3447  nn = sum(ncore(l, :)) + econfx(l)
3448  ii = 0
3449  DO
3450  ii = ii + 1
3451  IF (nn <= ll) THEN
3452  nelem(l, ii) = nn
3453  EXIT
3454  ELSE
3455  nelem(l, ii) = ll
3456  nn = nn - ll
3457  END IF
3458  END DO
3459  END DO
3460  ncalc = nelem - ncore
3461  ELSE
3462  ! test for compatibility of valence occupation and full atomic occupation
3463  IF (iounit > 0) THEN
3464  WRITE (iounit, "(/,A,A2)") "WARNING: Core states irregular for atom type ", ptable(z)%symbol
3465  WRITE (iounit, "(A,10I3)") "WARNING: Redefine ELEC_CONF in the KIND section"
3466  cpabort("Incompatible Atomic Occupations Detected")
3467  END IF
3468  END IF
3469  ELSE
3470  lmin = min(lmat, ubound(ptable(z)%e_conv, 1))
3471  ncore = 0
3472  ncalc = 0
3473  DO l = 0, lmin
3474  ll = 2*(2*l + 1)
3475  nn = abs(econfx(l))
3476  ii = 0
3477  DO
3478  ii = ii + 1
3479  IF (nn <= ll) THEN
3480  ncalc(l, ii) = -nn
3481  EXIT
3482  ELSE
3483  ncalc(l, ii) = -ll
3484  nn = nn - ll
3485  END IF
3486  END DO
3487  END DO
3488  nelem = ncalc
3489  END IF
3490 
3491  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  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  nlcc = .false.
3509 
3510  DO ikind = 1, SIZE(qs_kind_set)
3511  CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential, sgp_potential=sgp_potential)
3512  IF (ASSOCIATED(gth_potential)) THEN
3513  CALL get_potential(potential=gth_potential, nlcc_present=nlcc_present)
3514  nlcc = nlcc .OR. nlcc_present
3515  ELSEIF (ASSOCIATED(sgp_potential)) THEN
3516  CALL get_potential(potential=sgp_potential, has_nlcc=nlcc_present)
3517  nlcc = nlcc .OR. nlcc_present
3518  END IF
3519  END DO
3520 
3521  END FUNCTION has_nlcc
3522 
3523 ! **************************************************************************************************
3524 
3525 END MODULE qs_kind_types
static int imax(int x, int y)
Returns the larger of two given integer (missing from the C standard)
subroutine, public atom_sgp_release(sgp_pot)
...
Definition: atom_sgp.F:767
subroutine, public sgp_construction(sgp_pot, ecp_pot, upf_pot, orb_basis, error)
...
Definition: atom_sgp.F:79
Define the atom type and its sub types.
Definition: atom_types.F:15
integer, parameter, public lmat
Definition: atom_types.F:67
subroutine, public read_ecp_potential(element_symbol, potential, pseudo_name, pseudo_file, potential_section)
...
Definition: atom_types.F:2841
Routines that process Quantum Espresso UPF files.
Definition: atom_upf.F:14
subroutine, public atom_read_upf(pot, upf_filename, read_header)
...
Definition: atom_upf.F:102
pure subroutine, public atom_release_upf(upfpot)
...
Definition: atom_upf.F:875
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
subroutine, public remove_basis_set_container(basis)
...
subroutine, public remove_basis_from_container(container, inum, basis_type)
...
subroutine, public get_basis_from_container(container, basis_set, inumbas, basis_type)
Retrieve a basis set from the container.
subroutine, public add_basis_set_to_container(container, basis_set, basis_set_type)
...
subroutine, public write_orb_basis_set(orb_basis_set, output_unit, header)
Write a Gaussian-type orbital (GTO) basis set data set to the output unit.
subroutine, public deallocate_sto_basis_set(sto_basis_set)
...
subroutine, public init_aux_basis_set(gto_basis_set)
...
subroutine, public allocate_gto_basis_set(gto_basis_set)
...
subroutine, public combine_basis_sets(basis_set, basis_set_add)
...
subroutine, public write_gto_basis_set(gto_basis_set, output_unit, header)
Write a Gaussian-type orbital (GTO) basis set data set to the output unit.
subroutine, public allocate_sto_basis_set(sto_basis_set)
...
subroutine, public create_gto_from_sto_basis(sto_basis_set, gto_basis_set, ngauss, ortho)
...
subroutine, public read_sto_basis_set(element_symbol, basis_set_name, sto_basis_set, para_env, dft_section)
...
subroutine, public init_orb_basis_set(gto_basis_set)
Initialise a Gaussian-type orbital (GTO) basis set data set.
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius)
...
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
various routines to log and control the output. The idea is that decisions about where to log should ...
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
integer, parameter, public cp_p_file
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition of the atomic potential types.
subroutine, public set_default_all_potential(potential, z, zeff_correction)
...
subroutine, public create_1c_basis(orb_basis, soft_basis, gapw_1c_basis, basis_1c_level)
create the one center basis from the orbital basis
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public gapw_1c_large
integer, parameter, public do_method_pdg
integer, parameter, public do_method_pnnl
integer, parameter, public gapw_1c_medium
integer, parameter, public do_method_pw
integer, parameter, public do_method_rm1
integer, parameter, public gapw_1c_small
integer, parameter, public do_method_pm3
integer, parameter, public do_sirius
integer, parameter, public do_method_mndo
integer, parameter, public gapw_1c_orb
integer, parameter, public gapw_1c_very_large
integer, parameter, public do_method_mndod
integer, parameter, public do_method_am1
integer, parameter, public do_method_dftb
integer, parameter, public do_qs
integer, parameter, public do_method_xtb
integer, parameter, public do_method_pm6fm
integer, parameter, public do_method_pm6
objects that represent the structure of input sections and the data contained in an input section
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
integer, parameter, public default_string_length
Definition: kinds.F:57
integer, parameter, public default_path_length
Definition: kinds.F:58
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
real(kind=dp), parameter, public pi
Interface to the message passing library MPI.
Provides Cartesian and spherical orbital pointers and indices.
subroutine, public init_orbital_pointers(maxl)
Initialize or update the orbital pointers.
integer, dimension(:), allocatable, public nco
integer, dimension(:), allocatable, public ncoset
Factory routines for potentials used e.g. by pao_param_exp and pao_ml.
subroutine, public get_paw_proj_set(paw_proj_set, csprj, chprj, first_prj, first_prjs, last_prj, local_oce_sphi_h, local_oce_sphi_s, maxl, ncgauprj, nsgauprj, nsatbas, nsotot, nprj, o2nindex, n2oindex, rcprj, rzetprj, zisomin, zetprj)
Get informations about a paw projectors set.
subroutine, public allocate_paw_proj_set(paw_proj_set)
Allocate projector type for GAPW.
subroutine, public projectors(paw_proj, basis_1c, orb_basis, rc, qs_control, max_rad_local_type, force_env_section)
Initialize the projector-type set data set.
subroutine, public deallocate_paw_proj_set(paw_proj_set)
Deallocate a projector-type set data set.
Periodic Table related data definitions.
type(atom), dimension(0:nelem), public ptable
subroutine, public get_ptable_info(symbol, number, amass, ielement, covalent_radius, metallic_radius, vdw_radius, found)
Pass information about the kind given the element symbol.
Definition of physical constants:
Definition: physcon.F:68
real(kind=dp), parameter, public evolt
Definition: physcon.F:183
real(kind=dp), parameter, public angstrom
Definition: physcon.F:144
real(kind=dp), parameter, public bohr
Definition: physcon.F:147
Definition of the DFTB parameter types.
Definition: qs_dftb_types.F:12
Working with the DFTB parameter types.
Definition: qs_dftb_utils.F:12
subroutine, public deallocate_dftb_atom_param(dftb_parameter)
...
Definition: qs_dftb_utils.F:93
subroutine, public write_dftb_atom_param(dftb_parameter, subsys_section)
...
subroutine, public get_dftb_atom_param(dftb_parameter, name, typ, defined, z, zeff, natorb, lmax, skself, occupation, eta, energy, cutoff, xi, di, rcdisp, dudq)
...
Definition of disperson types for DFT calculations.
subroutine, public deallocate_grid_atom(grid_atom)
Deallocate a Gaussian-type orbital (GTO) basis set data set.
Definition: qs_grid_atom.F:105
subroutine, public allocate_grid_atom(grid_atom)
Initialize components of the grid_atom_type structure.
Definition: qs_grid_atom.F:73
subroutine, public allocate_harmonics_atom(harmonics)
Allocate a spherical harmonics set for the atom grid.
subroutine, public deallocate_harmonics_atom(harmonics)
Deallocate the spherical harmonics set for the atom grid.
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
subroutine, public set_pseudo_state(econf, z, ncalc, ncore, nelem)
...
logical function, public has_nlcc(qs_kind_set)
finds if a given qs run needs to use nlcc
subroutine, public deallocate_qs_kind_set(qs_kind_set)
Destructor routine for a set of qs kinds.
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
subroutine, public init_atom_electronic_state(atomic_kind, qs_kind, ncalc, ncore, nelem, edelta)
...
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr)
Get attributes of an atomic kind set.
subroutine, public init_gapw_nlcc(qs_kind_set)
...
subroutine, public set_qs_kind(qs_kind, paw_atom, ghost, floating, hard_radius, hard0_radius, covalent_radius, vdw_radius, lmax_rho0, zeff, no_optimize, dispersion, u_minus_j, reltmat, dftb_parameter, xtb_parameter, elec_conf, pao_basis_size)
Set the components of an atomic kind data set.
subroutine, public write_qs_kind_set(qs_kind_set, subsys_section)
Write an atomic kind set data set to the output unit.
subroutine, public init_gapw_basis_set(qs_kind_set, qs_control, force_env_section, modify_qs_control)
...
subroutine, public create_qs_kind_set(qs_kind_set, atomic_kind_set, kind_section, para_env, force_env_section)
Read an atomic kind set data set from the input file.
subroutine, public init_qs_kind_set(qs_kind_set)
Initialise an atomic kind set data set.
subroutine, public check_qs_kind_set(qs_kind_set, dft_control, subsys_section)
...
subroutine, public write_gto_basis_sets(qs_kind_set, subsys_section)
Write all the GTO basis sets of an atomic kind set to the output unit (for the printing of the unnorm...
Definition of the semi empirical parameter types.
subroutine, public write_se_param(sep, subsys_section)
Writes the semi-empirical type.
subroutine, public semi_empirical_create(sep)
Allocate semi-empirical type.
subroutine, public get_se_param(sep, name, typ, defined, z, zeff, natorb, eheat, beta, sto_exponents, uss, upp, udd, uff, alp, eisol, gss, gsp, gpp, gp2, acoul, nr, de, ass, asp, app, hsp, gsd, gpd, gdd, ppddg, dpddg, ngauss)
Get info from the semi-empirical type.
subroutine, public semi_empirical_release(sep)
Deallocate the semi-empirical type.
Working with the semi empirical parameter types.
subroutine, public se_param_set_default(sep, z, method)
Initialize parameter for a semi_empirival type.
subroutine, public init_se_param(sep, orb_basis_set, ngauss)
Initialize semi_empirical type.
subroutine, public create_soft_basis(orb_basis, soft_basis, eps_fit, rc, paw_atom, paw_type_forced, gpw_r3d_rs_type_forced)
create the soft basis from a GTO basis
Utilities for string manipulations.
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.
Read xTB parameters.
real(kind=dp) function, public xtb_set_kab(za, zb, xtb_control)
...
Definition of the xTB parameter types.
Definition: xtb_types.F:20
subroutine, public get_xtb_atom_param(xtb_parameter, symbol, aname, typ, defined, z, zeff, natorb, lmax, nao, lao, rcut, rcov, kx, eta, xgamma, alpha, zneff, nshell, nval, lval, kpoly, kappa, hen, zeta, occupation, electronegativity, chmax)
...
Definition: xtb_types.F:175
subroutine, public deallocate_xtb_atom_param(xtb_parameter)
...
Definition: xtb_types.F:133
subroutine, public write_xtb_atom_param(xtb_parameter, subsys_section)
...
Definition: xtb_types.F:313