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