(git:6344619)
Loading...
Searching...
No Matches
atom_types.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Define the atom type and its sub types
10!> \author jgh
11!> \date 03.03.2008
12!> \version 1.0
13!>
14! **************************************************************************************************
16 USE atom_upf, ONLY: atom_read_upf,&
18 atom_upfpot_type
19 USE bessel_lib, ONLY: bessel0
20 USE bibliography, ONLY: limpanuparb2011,&
21 cite_reference
32 USE input_constants, ONLY: &
42 USE input_val_types, ONLY: val_get,&
44 USE kinds, ONLY: default_path_length,&
46 dp
47 USE mathconstants, ONLY: dfac,&
48 fac,&
49 pi,&
50 rootpi
52 ptable
57 USE string_utilities, ONLY: remove_word,&
59#include "./base/base_uses.f90"
60
61 IMPLICIT NONE
62
63 PRIVATE
64
65 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_types'
66
67 ! maximum l-quantum number considered in atomic code/basis
68 INTEGER, PARAMETER :: lmat = 5
69
70 INTEGER, PARAMETER :: gto_basis = 100, &
71 cgto_basis = 101, &
72 sto_basis = 102, &
73 num_basis = 103
74
75 INTEGER, PARAMETER :: nmax = 25
76
77!> \brief Provides all information about a basis set
78! **************************************************************************************************
80 INTEGER :: basis_type = gto_basis
81 INTEGER, DIMENSION(0:lmat) :: nbas = 0
82 INTEGER, DIMENSION(0:lmat) :: nprim = 0
83 REAL(kind=dp), DIMENSION(:, :), POINTER :: am => null() !GTO exponents
84 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: cm => null() !Contraction coeffs
85 REAL(kind=dp), DIMENSION(:, :), POINTER :: as => null() !STO exponents
86 INTEGER, DIMENSION(:, :), POINTER :: ns => null() !STO n-quantum numbers
87 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: bf => null() !num. bsf
88 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: dbf => null() !derivatives (num)
89 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: ddbf => null() !2nd derivatives (num)
90 REAL(kind=dp) :: eps_eig = 0.0_dp
91 TYPE(grid_atom_type), POINTER :: grid => null()
92 LOGICAL :: geometrical = .false.
93 REAL(kind=dp) :: aval = 0.0_dp, cval = 0.0_dp
94 INTEGER, DIMENSION(0:lmat) :: start = 0
95 END TYPE atom_basis_type
96
97!> \brief Provides all information about a pseudopotential
98! **************************************************************************************************
100 CHARACTER(LEN=2) :: symbol = ""
101 CHARACTER(LEN=default_string_length) :: pname = ""
102 INTEGER, DIMENSION(0:lmat) :: econf = 0
103 REAL(dp) :: zion = 0.0_dp
104 REAL(dp) :: rc = 0.0_dp
105 INTEGER :: ncl = 0
106 REAL(dp), DIMENSION(5) :: cl = 0.0_dp
107 INTEGER, DIMENSION(0:lmat) :: nl = 0
108 REAL(dp), DIMENSION(0:lmat) :: rcnl = 0.0_dp
109 REAL(dp), DIMENSION(4, 4, 0:lmat) :: hnl = 0.0_dp
110 ! SOC
111 LOGICAL :: soc = .false.
112 REAL(dp), DIMENSION(4, 4, 0:lmat) :: knl = 0.0_dp
113 ! type extensions
114 ! NLCC
115 LOGICAL :: nlcc = .false.
116 INTEGER :: nexp_nlcc = 0
117 REAL(kind=dp), DIMENSION(10) :: alpha_nlcc = 0.0_dp
118 INTEGER, DIMENSION(10) :: nct_nlcc = 0
119 REAL(kind=dp), DIMENSION(4, 10) :: cval_nlcc = 0.0_dp
120 ! LSD potential
121 LOGICAL :: lsdpot = .false.
122 INTEGER :: nexp_lsd = 0
123 REAL(kind=dp), DIMENSION(10) :: alpha_lsd = 0.0_dp
124 INTEGER, DIMENSION(10) :: nct_lsd = 0
125 REAL(kind=dp), DIMENSION(4, 10) :: cval_lsd = 0.0_dp
126 ! extended local potential
127 LOGICAL :: lpotextended = .false.
128 INTEGER :: nexp_lpot = 0
129 REAL(kind=dp), DIMENSION(10) :: alpha_lpot = 0.0_dp
130 INTEGER, DIMENSION(10) :: nct_lpot = 0
131 REAL(kind=dp), DIMENSION(4, 10) :: cval_lpot = 0.0_dp
132 END TYPE atom_gthpot_type
133
135 CHARACTER(LEN=2) :: symbol = ""
136 CHARACTER(LEN=default_string_length) :: pname = ""
137 INTEGER, DIMENSION(0:lmat) :: econf = 0
138 REAL(dp) :: zion = 0.0_dp
139 INTEGER :: lmax = 0
140 INTEGER :: nloc = 0 ! # terms
141 INTEGER, DIMENSION(1:15) :: nrloc = 0 ! r**(n-2)
142 REAL(dp), DIMENSION(1:15) :: aloc = 0.0_dp ! coefficient
143 REAL(dp), DIMENSION(1:15) :: bloc = 0.0_dp ! exponent
144 INTEGER, DIMENSION(0:10) :: npot = 0 ! # terms
145 INTEGER, DIMENSION(1:15, 0:10) :: nrpot = 0 ! r**(n-2)
146 REAL(dp), DIMENSION(1:15, 0:10) :: apot = 0.0_dp ! coefficient
147 REAL(dp), DIMENSION(1:15, 0:10) :: bpot = 0.0_dp ! exponent
148 END TYPE atom_ecppot_type
149
151 CHARACTER(LEN=2) :: symbol = ""
152 CHARACTER(LEN=default_string_length) :: pname = ""
153 INTEGER, DIMENSION(0:lmat) :: econf = 0
154 REAL(dp) :: zion = 0.0_dp
155 INTEGER :: lmax = 0
156 LOGICAL :: has_nonlocal = .false.
157 INTEGER :: n_nonlocal = 0
158 LOGICAL, DIMENSION(0:5) :: is_nonlocal = .false.
159 REAL(kind=dp), DIMENSION(nmax) :: a_nonlocal = 0.0_dp
160 REAL(kind=dp), DIMENSION(nmax, 0:lmat) :: h_nonlocal = 0.0_dp
161 REAL(kind=dp), DIMENSION(nmax, nmax, 0:lmat) :: c_nonlocal = 0.0_dp
162 INTEGER :: n_local = 0
163 REAL(kind=dp) :: ac_local = 0.0_dp
164 REAL(kind=dp), DIMENSION(nmax) :: a_local = 0.0_dp
165 REAL(kind=dp), DIMENSION(nmax) :: c_local = 0.0_dp
166 LOGICAL :: has_nlcc = .false.
167 INTEGER :: n_nlcc = 0
168 REAL(kind=dp), DIMENSION(nmax) :: a_nlcc = 0.0_dp
169 REAL(kind=dp), DIMENSION(nmax) :: c_nlcc = 0.0_dp
170 END TYPE atom_sgppot_type
171
173 INTEGER :: ppot_type = 0
174 LOGICAL :: confinement = .false.
175 INTEGER :: conf_type = 0
176 REAL(dp) :: acon = 0.0_dp
177 REAL(dp) :: rcon = 0.0_dp
178 REAL(dp) :: scon = 0.0_dp
181 TYPE(atom_upfpot_type) :: upf_pot = atom_upfpot_type()
183 END TYPE atom_potential_type
184
185!> \brief Provides info about hartree-fock exchange (For now, we only support potentials that can be represented
186!> with Coulomb and longrange-coulomb potential)
187! **************************************************************************************************
189 REAL(kind=dp) :: scale_coulomb = 0.0_dp
190 REAL(kind=dp) :: scale_longrange = 0.0_dp
191 REAL(kind=dp) :: omega = 0.0_dp
192 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: kernel
193 LOGICAL :: do_gh = .false.
194 INTEGER :: nr_gh = 0
195 END TYPE atom_hfx_type
196
197!> \brief Provides all information on states and occupation
198! **************************************************************************************************
200 REAL(kind=dp), DIMENSION(0:lmat, 10) :: occ = 0.0_dp
201 REAL(kind=dp), DIMENSION(0:lmat, 10) :: core = 0.0_dp
202 REAL(kind=dp), DIMENSION(0:lmat, 10) :: occupation = 0.0_dp
203 INTEGER :: maxl_occ = 0
204 INTEGER, DIMENSION(0:lmat) :: maxn_occ = 0
205 INTEGER :: maxl_calc = 0
206 INTEGER, DIMENSION(0:lmat) :: maxn_calc = 0
207 INTEGER :: multiplicity = 0
208 REAL(kind=dp), DIMENSION(0:lmat, 10) :: occa = 0.0_dp, occb = 0.0_dp
209 END TYPE atom_state
210
211!> \brief Holds atomic integrals
212! **************************************************************************************************
213 TYPE eri
214 REAL(kind=dp), DIMENSION(:, :), POINTER :: int => null()
215 END TYPE eri
216
218 INTEGER :: status = 0
219 INTEGER :: ppstat = 0
220 LOGICAL :: eri_coulomb = .false.
221 LOGICAL :: eri_exchange = .false.
222 LOGICAL :: all_nu = .false.
223 INTEGER, DIMENSION(0:lmat) :: n = 0, nne = 0
224 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: ovlp => null(), kin => null(), core => null(), clsd => null()
225 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: utrans => null(), uptrans => null()
226 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: hnl => null()
227 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: conf => null()
228 TYPE(eri), DIMENSION(100) :: ceri = eri()
229 TYPE(eri), DIMENSION(100) :: eeri = eri()
230 INTEGER :: dkhstat = 0
231 INTEGER :: zorastat = 0
232 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: tzora => null()
233 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: hdkh => null()
234 END TYPE atom_integrals
235
236!> \brief Holds atomic orbitals and energies
237! **************************************************************************************************
239 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: wfn => null(), wfna => null(), wfnb => null()
240 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: pmat => null(), pmata => null(), pmatb => null()
241 REAL(kind=dp), DIMENSION(:, :), POINTER :: ener => null(), enera => null(), enerb => null()
242 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: refene => null(), refchg => null(), refnod => null()
243 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: wrefene => null(), wrefchg => null(), wrefnod => null()
244 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: crefene => null(), crefchg => null(), crefnod => null()
245 REAL(kind=dp), DIMENSION(:, :), POINTER :: wpsir0 => null(), tpsir0 => null()
246 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: rcmax => null()
247 CHARACTER(LEN=2), DIMENSION(:, :, :), POINTER :: reftype => null()
248 END TYPE atom_orbitals
249
250!> \brief Operator matrices
251! **************************************************************************************************
253 INTEGER, DIMENSION(0:lmat) :: n = 0
254 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: op => null()
255 END TYPE opmat_type
256
257!> \brief Operator grids
258! **************************************************************************************************
260 REAL(kind=dp), DIMENSION(:), POINTER :: op => null()
261 TYPE(grid_atom_type), POINTER :: grid => null()
262 END TYPE opgrid_type
263
264!> \brief All energies
265! **************************************************************************************************
266 TYPE atom_energy_type
267 REAL(kind=dp) :: etot = 0.0_dp
268 REAL(kind=dp) :: eband = 0.0_dp
269 REAL(kind=dp) :: ekin = 0.0_dp
270 REAL(kind=dp) :: epot = 0.0_dp
271 REAL(kind=dp) :: ecore = 0.0_dp
272 REAL(kind=dp) :: elsd = 0.0_dp
273 REAL(kind=dp) :: epseudo = 0.0_dp
274 REAL(kind=dp) :: eploc = 0.0_dp
275 REAL(kind=dp) :: epnl = 0.0_dp
276 REAL(kind=dp) :: exc = 0.0_dp
277 REAL(kind=dp) :: ecoulomb = 0.0_dp
278 REAL(kind=dp) :: eexchange = 0.0_dp
279 REAL(kind=dp) :: econfinement = 0.0_dp
280 END TYPE atom_energy_type
281
282!> \brief Information on optimization procedure
283! **************************************************************************************************
285 REAL(kind=dp) :: damping = 0.0_dp
286 REAL(kind=dp) :: eps_scf = 0.0_dp
287 REAL(kind=dp) :: eps_diis = 0.0_dp
288 INTEGER :: max_iter = 0
289 INTEGER :: n_diis = 0
291
292!> \brief Provides all information about an atomic kind
293! **************************************************************************************************
295 INTEGER :: z = 0
296 INTEGER :: zcore = 0
297 LOGICAL :: pp_calc = .false.
298! ZMP adding in type some variables
299 LOGICAL :: do_zmp = .false., doread = .false., read_vxc = .false., dm = .false.
300 CHARACTER(LEN=default_string_length) :: ext_file = "", ext_vxc_file = "", &
301 zmp_restart_file = ""
302!
303 INTEGER :: method_type = do_rks_atom
304 INTEGER :: relativistic = do_nonrel_atom
305 INTEGER :: coulomb_integral_type = do_analytic
306 INTEGER :: exchange_integral_type = do_analytic
307! ZMP
308 REAL(kind=dp) :: lambda = 0.0_dp
309 REAL(kind=dp) :: rho_diff_integral = 0.0_dp
310 REAL(kind=dp) :: weight = 0.0_dp, zmpgrid_tol = 0.0_dp, zmpvxcgrid_tol = 0.0_dp
311!
312 TYPE(atom_basis_type), POINTER :: basis => null()
313 TYPE(atom_potential_type), POINTER :: potential => null()
314 TYPE(atom_state), POINTER :: state => null()
315 TYPE(atom_integrals), POINTER :: integrals => null()
316 TYPE(atom_orbitals), POINTER :: orbitals => null()
317 TYPE(atom_energy_type) :: energy = atom_energy_type()
319 TYPE(section_vals_type), POINTER :: xc_section => null(), zmp_section => null()
320 TYPE(opmat_type), POINTER :: fmat => null()
321 TYPE(atom_hfx_type) :: hfx_pot = atom_hfx_type()
322 END TYPE atom_type
323! **************************************************************************************************
325 TYPE(atom_type), POINTER :: atom => null()
326 END TYPE atom_p_type
327
328 PUBLIC :: lmat
332 PUBLIC :: atom_optimization_type
333 PUBLIC :: atom_compare_grids
339 PUBLIC :: clementi_geobas
344 PUBLIC :: setup_hf_section
345
347 MODULE PROCEDURE read_ecp_potential_file, &
348 read_ecp_potential_files
349 END INTERFACE
350
351! **************************************************************************************************
352
353CONTAINS
354
355! **************************************************************************************************
356!> \brief Initialize the basis for the atomic code
357!> \param basis ...
358!> \param basis_section ...
359!> \param zval ...
360!> \param btyp ...
361!> \note Highly accurate relativistic universal Gaussian basis set: Dirac-Fock-Coulomb calculations
362!> for atomic systems up to nobelium
363!> J. Chem. Phys. 101, 6829 (1994); DOI:10.1063/1.468311
364!> G. L. Malli and A. B. F. Da Silva
365!> Department of Chemistry, Simon Fraser University, Burnaby, B.C., Canada
366!> Yasuyuki Ishikawa
367!> Department of Chemistry, University of Puerto Rico, San Juan, Puerto Rico
368!>
369!> A universal Gaussian basis set is developed that leads to relativistic Dirac-Fock SCF energies
370!> of comparable accuracy as that obtained by the accurate numerical finite-difference method
371!> (GRASP2 package) [J. Phys. B 25, 1 (1992)]. The Gaussian-type functions of our universal basis
372!> set satisfy the relativistic boundary conditions associated with the finite nuclear model for a
373!> finite speed of light and conform to the so-called kinetic balance at the nonrelativistic limit.
374!> We attribute the exceptionally high accuracy obtained in our calculations to the fact that the
375!> representation of the relativistic dynamics of an electron in a spherical ball finite nucleus
376!> near the origin in terms of our universal Gaussian basis set is as accurate as that provided by
377!> the numerical finite-difference method. Results of the Dirac-Fock-Coulomb energies for a number
378!> of atoms up to No (Z=102) and some negative ions are presented and compared with the recent
379!> results obtained with the numerical finite-difference method and geometrical Gaussian basis sets
380!> by Parpia, Mohanty, and Clementi [J. Phys. B 25, 1 (1992)]. The accuracy of our calculations is
381!> estimated to be within a few parts in 109 for all the atomic systems studied.
382! **************************************************************************************************
383 SUBROUTINE init_atom_basis(basis, basis_section, zval, btyp)
384 TYPE(atom_basis_type), INTENT(INOUT) :: basis
385 TYPE(section_vals_type), POINTER :: basis_section
386 INTEGER, INTENT(IN) :: zval
387 CHARACTER(LEN=2) :: btyp
388
389 CHARACTER(LEN=*), PARAMETER :: routinen = 'init_atom_basis'
390 INTEGER, PARAMETER :: nua = 40, nup = 16
391 REAL(kind=dp), DIMENSION(nua), PARAMETER :: ugbs = [0.007299_dp, 0.013705_dp, 0.025733_dp, &
392 0.048316_dp, 0.090718_dp, 0.170333_dp, 0.319819_dp, 0.600496_dp, 1.127497_dp, 2.117000_dp,&
393 3.974902_dp, 7.463317_dp, 14.013204_dp, 26.311339_dp, 49.402449_dp, 92.758561_dp, &
394 174.164456_dp, 327.013024_dp, 614.003114_dp, 1152.858743_dp, 2164.619772_dp, &
395 4064.312984_dp, 7631.197056_dp, 14328.416324_dp, 26903.186074_dp, 50513.706789_dp, &
396 94845.070265_dp, 178082.107320_dp, 334368.848683_dp, 627814.487663_dp, 1178791.123851_dp, &
397 2213310.684886_dp, 4155735.557141_dp, 7802853.046713_dp, 14650719.428954_dp, &
398 27508345.793637_dp, 51649961.080194_dp, 96978513.342764_dp, 182087882.613702_dp, &
399 341890134.751331_dp]
400
401 CHARACTER(LEN=default_path_length) :: basis_fn
402 CHARACTER(LEN=default_string_length) :: basis_name
403 INTEGER :: basistype, handle, i, j, k, l, ll, m, &
404 ngp, nl, nr, nu, quadtype
405 INTEGER, DIMENSION(0:lmat) :: starti
406 INTEGER, DIMENSION(:), POINTER :: nqm, num_gto, num_slater, sindex
407 REAL(kind=dp) :: al, amax, aval, cval, ear, pf, rk
408 REAL(kind=dp), DIMENSION(:), POINTER :: expo
409 TYPE(section_vals_type), POINTER :: gto_basis_section
410
411 CALL timeset(routinen, handle)
412
413 ! btyp = AE : standard all-electron basis
414 ! btyp = PP : standard pseudopotential basis
415 ! btyp = AA : high accuracy all-electron basis
416 ! btyp = AP : high accuracy pseudopotential basis
417
418 NULLIFY (basis%am, basis%cm, basis%as, basis%ns, basis%bf, basis%dbf, basis%ddbf)
419 ! get information on quadrature type and number of grid points
420 ! allocate and initialize the atomic grid
421 CALL allocate_grid_atom(basis%grid)
422 CALL section_vals_val_get(basis_section, "QUADRATURE", i_val=quadtype)
423 CALL section_vals_val_get(basis_section, "GRID_POINTS", i_val=ngp)
424 IF (ngp <= 0) THEN
425 cpabort("The number of radial grid points must be greater than zero.")
426 END IF
427 CALL create_grid_atom(basis%grid, ngp, 1, 1, 0, quadtype)
428 basis%grid%nr = ngp
429 basis%geometrical = .false.
430 basis%aval = 0._dp
431 basis%cval = 0._dp
432 basis%start = 0
433
434 CALL section_vals_val_get(basis_section, "BASIS_TYPE", i_val=basistype)
435 CALL section_vals_val_get(basis_section, "EPS_EIGENVALUE", r_val=basis%eps_eig)
436 SELECT CASE (basistype)
437 CASE (gaussian)
438 basis%basis_type = gto_basis
439 NULLIFY (num_gto)
440 CALL section_vals_val_get(basis_section, "NUM_GTO", i_vals=num_gto)
441 IF (num_gto(1) < 1) THEN
442 ! use default basis
443 IF (btyp == "AE") THEN
444 nu = nua
445 ELSE IF (btyp == "PP") THEN
446 nu = nup
447 ELSE
448 nu = nua
449 END IF
450 basis%nbas = nu
451 basis%nprim = nu
452 ALLOCATE (basis%am(nu, 0:lmat))
453 DO i = 0, lmat
454 basis%am(1:nu, i) = ugbs(1:nu)
455 END DO
456 ELSE
457 basis%nbas = 0
458 DO i = 1, SIZE(num_gto)
459 basis%nbas(i - 1) = num_gto(i)
460 END DO
461 basis%nprim = basis%nbas
462 m = maxval(basis%nbas)
463 ALLOCATE (basis%am(m, 0:lmat))
464 basis%am = 0._dp
465 DO l = 0, lmat
466 IF (basis%nbas(l) > 0) THEN
467 NULLIFY (expo)
468 SELECT CASE (l)
469 CASE (0)
470 CALL section_vals_val_get(basis_section, "S_EXPONENTS", r_vals=expo)
471 CASE (1)
472 CALL section_vals_val_get(basis_section, "P_EXPONENTS", r_vals=expo)
473 CASE (2)
474 CALL section_vals_val_get(basis_section, "D_EXPONENTS", r_vals=expo)
475 CASE (3)
476 CALL section_vals_val_get(basis_section, "F_EXPONENTS", r_vals=expo)
477 CASE DEFAULT
478 cpabort("Invalid angular quantum number l found for Gaussian basis set")
479 END SELECT
480 cpassert(SIZE(expo) >= basis%nbas(l))
481 DO i = 1, basis%nbas(l)
482 basis%am(i, l) = expo(i)
483 END DO
484 END IF
485 END DO
486 END IF
487 ! initialize basis function on a radial grid
488 nr = basis%grid%nr
489 m = maxval(basis%nbas)
490 ALLOCATE (basis%bf(nr, m, 0:lmat))
491 ALLOCATE (basis%dbf(nr, m, 0:lmat))
492 ALLOCATE (basis%ddbf(nr, m, 0:lmat))
493 basis%bf = 0._dp
494 basis%dbf = 0._dp
495 basis%ddbf = 0._dp
496 DO l = 0, lmat
497 DO i = 1, basis%nbas(l)
498 al = basis%am(i, l)
499 DO k = 1, nr
500 rk = basis%grid%rad(k)
501 ear = exp(-al*basis%grid%rad(k)**2)
502 basis%bf(k, i, l) = rk**l*ear
503 basis%dbf(k, i, l) = (real(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear
504 basis%ddbf(k, i, l) = (real(l*(l - 1), dp)*rk**(l - 2) - &
505 2._dp*al*real(2*l + 1, dp)*rk**(l) + 4._dp*al**2*rk**(l + 2))*ear
506 END DO
507 END DO
508 END DO
509 CASE (geometrical_gto)
510 basis%basis_type = gto_basis
511 NULLIFY (num_gto)
512 CALL section_vals_val_get(basis_section, "NUM_GTO", i_vals=num_gto)
513 IF (num_gto(1) < 1) THEN
514 IF (btyp == "AE") THEN
515 ! use the Clementi extra large basis
516 CALL clementi_geobas(zval, cval, aval, basis%nbas, starti)
517 ELSE IF (btyp == "PP") THEN
518 ! use the Clementi extra large basis
519 CALL clementi_geobas(zval, cval, aval, basis%nbas, starti)
520 ELSE IF (btyp == "AA") THEN
521 CALL clementi_geobas(zval, cval, aval, basis%nbas, starti)
522 amax = cval**(basis%nbas(0) - 1)
523 basis%nbas(0) = nint((log(amax)/log(1.6_dp)))
524 cval = 1.6_dp
525 starti = 0
526 basis%nbas(1) = basis%nbas(0) - 4
527 basis%nbas(2) = basis%nbas(0) - 8
528 basis%nbas(3) = basis%nbas(0) - 12
529 IF (lmat > 3) basis%nbas(4:lmat) = 0
530 ELSE IF (btyp == "AP") THEN
531 CALL clementi_geobas(zval, cval, aval, basis%nbas, starti)
532 amax = 500._dp/aval
533 basis%nbas = nint((log(amax)/log(1.6_dp)))
534 cval = 1.6_dp
535 starti = 0
536 ELSE
537 ! use the Clementi extra large basis
538 CALL clementi_geobas(zval, cval, aval, basis%nbas, starti)
539 END IF
540 basis%nprim = basis%nbas
541 ELSE
542 basis%nbas = 0
543 DO i = 1, SIZE(num_gto)
544 basis%nbas(i - 1) = num_gto(i)
545 END DO
546 basis%nprim = basis%nbas
547 NULLIFY (sindex)
548 CALL section_vals_val_get(basis_section, "START_INDEX", i_vals=sindex)
549 starti = 0
550 DO i = 1, SIZE(sindex)
551 starti(i - 1) = sindex(i)
552 cpassert(sindex(i) >= 0)
553 END DO
554 CALL section_vals_val_get(basis_section, "GEOMETRICAL_FACTOR", r_val=cval)
555 CALL section_vals_val_get(basis_section, "GEO_START_VALUE", r_val=aval)
556 END IF
557 m = maxval(basis%nbas)
558 ALLOCATE (basis%am(m, 0:lmat))
559 basis%am = 0._dp
560 DO l = 0, lmat
561 DO i = 1, basis%nbas(l)
562 ll = i - 1 + starti(l)
563 basis%am(i, l) = aval*cval**(ll)
564 END DO
565 END DO
566
567 basis%geometrical = .true.
568 basis%aval = aval
569 basis%cval = cval
570 basis%start = starti
571
572 ! initialize basis function on a radial grid
573 nr = basis%grid%nr
574 m = maxval(basis%nbas)
575 ALLOCATE (basis%bf(nr, m, 0:lmat))
576 ALLOCATE (basis%dbf(nr, m, 0:lmat))
577 ALLOCATE (basis%ddbf(nr, m, 0:lmat))
578 basis%bf = 0._dp
579 basis%dbf = 0._dp
580 basis%ddbf = 0._dp
581 DO l = 0, lmat
582 DO i = 1, basis%nbas(l)
583 al = basis%am(i, l)
584 DO k = 1, nr
585 rk = basis%grid%rad(k)
586 ear = exp(-al*basis%grid%rad(k)**2)
587 basis%bf(k, i, l) = rk**l*ear
588 basis%dbf(k, i, l) = (real(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear
589 basis%ddbf(k, i, l) = (real(l*(l - 1), dp)*rk**(l - 2) - &
590 2._dp*al*real(2*l + 1, dp)*rk**(l) + 4._dp*al**2*rk**(l + 2))*ear
591 END DO
592 END DO
593 END DO
594 CASE (contracted_gto)
595 basis%basis_type = cgto_basis
596 CALL section_vals_val_get(basis_section, "BASIS_SET_FILE_NAME", c_val=basis_fn)
597 CALL section_vals_val_get(basis_section, "BASIS_SET", c_val=basis_name)
598 gto_basis_section => section_vals_get_subs_vals(basis_section, "BASIS")
599 CALL read_basis_set(ptable(zval)%symbol, basis, basis_name, basis_fn, &
600 gto_basis_section)
601
602 ! initialize basis function on a radial grid
603 nr = basis%grid%nr
604 m = maxval(basis%nbas)
605 ALLOCATE (basis%bf(nr, m, 0:lmat))
606 ALLOCATE (basis%dbf(nr, m, 0:lmat))
607 ALLOCATE (basis%ddbf(nr, m, 0:lmat))
608 basis%bf = 0._dp
609 basis%dbf = 0._dp
610 basis%ddbf = 0._dp
611 DO l = 0, lmat
612 DO i = 1, basis%nprim(l)
613 al = basis%am(i, l)
614 DO k = 1, nr
615 rk = basis%grid%rad(k)
616 ear = exp(-al*basis%grid%rad(k)**2)
617 DO j = 1, basis%nbas(l)
618 basis%bf(k, j, l) = basis%bf(k, j, l) + rk**l*ear*basis%cm(i, j, l)
619 basis%dbf(k, j, l) = basis%dbf(k, j, l) &
620 + (real(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear*basis%cm(i, j, l)
621 basis%ddbf(k, j, l) = basis%ddbf(k, j, l) + &
622 (real(l*(l - 1), dp)*rk**(l - 2) - 2._dp*al*real(2*l + 1, dp)*rk**(l) + 4._dp*al**2*rk**(l + 2))* &
623 ear*basis%cm(i, j, l)
624 END DO
625 END DO
626 END DO
627 END DO
628 CASE (slater)
629 basis%basis_type = sto_basis
630 NULLIFY (num_slater)
631 CALL section_vals_val_get(basis_section, "NUM_SLATER", i_vals=num_slater)
632 IF (num_slater(1) < 1) THEN
633 cpabort("Invalid number (less than 1) Slater-type functions found.")
634 ELSE
635 basis%nbas = 0
636 DO i = 1, SIZE(num_slater)
637 basis%nbas(i - 1) = num_slater(i)
638 END DO
639 basis%nprim = basis%nbas
640 m = maxval(basis%nbas)
641 ALLOCATE (basis%as(m, 0:lmat), basis%ns(m, 0:lmat))
642 basis%as = 0.0_dp
643 basis%ns = 0
644 DO l = 0, lmat
645 IF (basis%nbas(l) > 0) THEN
646 NULLIFY (expo)
647 SELECT CASE (l)
648 CASE (0)
649 CALL section_vals_val_get(basis_section, "S_EXPONENTS", r_vals=expo)
650 CASE (1)
651 CALL section_vals_val_get(basis_section, "P_EXPONENTS", r_vals=expo)
652 CASE (2)
653 CALL section_vals_val_get(basis_section, "D_EXPONENTS", r_vals=expo)
654 CASE (3)
655 CALL section_vals_val_get(basis_section, "F_EXPONENTS", r_vals=expo)
656 CASE DEFAULT
657 cpabort("Invalid angular quantum number l found for Slater basis set")
658 END SELECT
659 cpassert(SIZE(expo) >= basis%nbas(l))
660 DO i = 1, basis%nbas(l)
661 basis%as(i, l) = expo(i)
662 END DO
663 NULLIFY (nqm)
664 SELECT CASE (l)
665 CASE (0)
666 CALL section_vals_val_get(basis_section, "S_QUANTUM_NUMBERS", i_vals=nqm)
667 CASE (1)
668 CALL section_vals_val_get(basis_section, "P_QUANTUM_NUMBERS", i_vals=nqm)
669 CASE (2)
670 CALL section_vals_val_get(basis_section, "D_QUANTUM_NUMBERS", i_vals=nqm)
671 CASE (3)
672 CALL section_vals_val_get(basis_section, "F_QUANTUM_NUMBERS", i_vals=nqm)
673 CASE DEFAULT
674 cpabort("Invalid angular quantum number l found for Slater basis set")
675 END SELECT
676 cpassert(SIZE(nqm) >= basis%nbas(l))
677 DO i = 1, basis%nbas(l)
678 basis%ns(i, l) = nqm(i)
679 END DO
680 END IF
681 END DO
682 END IF
683 ! initialize basis function on a radial grid
684 nr = basis%grid%nr
685 m = maxval(basis%nbas)
686 ALLOCATE (basis%bf(nr, m, 0:lmat))
687 ALLOCATE (basis%dbf(nr, m, 0:lmat))
688 ALLOCATE (basis%ddbf(nr, m, 0:lmat))
689 basis%bf = 0._dp
690 basis%dbf = 0._dp
691 basis%ddbf = 0._dp
692 DO l = 0, lmat
693 DO i = 1, basis%nbas(l)
694 al = basis%as(i, l)
695 nl = basis%ns(i, l)
696 pf = (2._dp*al)**nl*sqrt(2._dp*al/fac(2*nl))
697 DO k = 1, nr
698 rk = basis%grid%rad(k)
699 ear = rk**(nl - 1)*exp(-al*rk)
700 basis%bf(k, i, l) = pf*ear
701 basis%dbf(k, i, l) = pf*(real(nl - 1, dp)/rk - al)*ear
702 basis%ddbf(k, i, l) = pf*(real((nl - 2)*(nl - 1), dp)/rk/rk &
703 - al*real(2*(nl - 1), dp)/rk + al*al)*ear
704 END DO
705 END DO
706 END DO
707 CASE (numerical)
708 basis%basis_type = num_basis
709 cpabort("Numerical basis set type not yet implemented.")
710 CASE DEFAULT
711 cpabort("Unknown basis set type specified. Check the code!")
712 END SELECT
713
714 CALL timestop(handle)
715
716 END SUBROUTINE init_atom_basis
717
718! **************************************************************************************************
719!> \brief ...
720!> \param basis ...
721! **************************************************************************************************
723 TYPE(atom_basis_type), INTENT(INOUT) :: basis
724
725 CHARACTER(LEN=*), PARAMETER :: routinen = 'init_atom_basis_default_pp'
726 INTEGER, PARAMETER :: nua = 40, nup = 20
727 REAL(kind=dp), DIMENSION(nua), PARAMETER :: ugbs = [0.007299_dp, 0.013705_dp, 0.025733_dp, &
728 0.048316_dp, 0.090718_dp, 0.170333_dp, 0.319819_dp, 0.600496_dp, 1.127497_dp, 2.117000_dp,&
729 3.974902_dp, 7.463317_dp, 14.013204_dp, 26.311339_dp, 49.402449_dp, 92.758561_dp, &
730 174.164456_dp, 327.013024_dp, 614.003114_dp, 1152.858743_dp, 2164.619772_dp, &
731 4064.312984_dp, 7631.197056_dp, 14328.416324_dp, 26903.186074_dp, 50513.706789_dp, &
732 94845.070265_dp, 178082.107320_dp, 334368.848683_dp, 627814.487663_dp, 1178791.123851_dp, &
733 2213310.684886_dp, 4155735.557141_dp, 7802853.046713_dp, 14650719.428954_dp, &
734 27508345.793637_dp, 51649961.080194_dp, 96978513.342764_dp, 182087882.613702_dp, &
735 341890134.751331_dp]
736
737 INTEGER :: handle, i, k, l, m, ngp, nr, nu, quadtype
738 REAL(kind=dp) :: al, ear, rk
739
740 CALL timeset(routinen, handle)
741
742 NULLIFY (basis%am, basis%cm, basis%as, basis%ns, basis%bf, basis%dbf, basis%ddbf)
743
744 ! Allocate and initialize the atomic grid
745 NULLIFY (basis%grid)
746 CALL allocate_grid_atom(basis%grid)
747 quadtype = do_gapw_log
748 ngp = 500
749 CALL create_grid_atom(basis%grid, ngp, 1, 1, 0, quadtype)
750 basis%grid%nr = ngp
751 basis%geometrical = .false.
752 basis%aval = 0._dp
753 basis%cval = 0._dp
754 basis%start = 0
755 basis%eps_eig = 1.e-12_dp
756
757 basis%basis_type = gto_basis
758 nu = nup
759 basis%nbas = nu
760 basis%nprim = nu
761 ALLOCATE (basis%am(nu, 0:lmat))
762 DO i = 0, lmat
763 basis%am(1:nu, i) = ugbs(1:nu)
764 END DO
765 ! initialize basis function on a radial grid
766 nr = basis%grid%nr
767 m = maxval(basis%nbas)
768 ALLOCATE (basis%bf(nr, m, 0:lmat))
769 ALLOCATE (basis%dbf(nr, m, 0:lmat))
770 ALLOCATE (basis%ddbf(nr, m, 0:lmat))
771 basis%bf = 0._dp
772 basis%dbf = 0._dp
773 basis%ddbf = 0._dp
774 DO l = 0, lmat
775 DO i = 1, basis%nbas(l)
776 al = basis%am(i, l)
777 DO k = 1, nr
778 rk = basis%grid%rad(k)
779 ear = exp(-al*basis%grid%rad(k)**2)
780 basis%bf(k, i, l) = rk**l*ear
781 basis%dbf(k, i, l) = (real(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear
782 basis%ddbf(k, i, l) = (real(l*(l - 1), dp)*rk**(l - 2) - &
783 2._dp*al*real(2*l + 1, dp)*rk**(l) + 4._dp*al**2*rk**(l + 2))*ear
784 END DO
785 END DO
786 END DO
787
788 CALL timestop(handle)
789
790 END SUBROUTINE init_atom_basis_default_pp
791
792! **************************************************************************************************
793!> \brief ...
794!> \param basis ...
795!> \param gbasis ...
796!> \param r ...
797!> \param rab ...
798! **************************************************************************************************
799 SUBROUTINE atom_basis_gridrep(basis, gbasis, r, rab)
800 TYPE(atom_basis_type), INTENT(IN) :: basis
801 TYPE(atom_basis_type), INTENT(INOUT) :: gbasis
802 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: r, rab
803
804 INTEGER :: i, j, k, l, m, n1, n2, n3, ngp, nl, nr
805 REAL(kind=dp) :: al, ear, pf, rk
806
807 NULLIFY (gbasis%am, gbasis%cm, gbasis%as, gbasis%ns, gbasis%bf, gbasis%dbf, gbasis%ddbf)
808
809 ! copy basis info
810 gbasis%basis_type = basis%basis_type
811 gbasis%nbas(0:lmat) = basis%nbas(0:lmat)
812 gbasis%nprim(0:lmat) = basis%nprim(0:lmat)
813 IF (ASSOCIATED(basis%am)) THEN
814 n1 = SIZE(basis%am, 1)
815 n2 = SIZE(basis%am, 2)
816 ALLOCATE (gbasis%am(n1, 0:n2 - 1))
817 gbasis%am = basis%am
818 END IF
819 IF (ASSOCIATED(basis%cm)) THEN
820 n1 = SIZE(basis%cm, 1)
821 n2 = SIZE(basis%cm, 2)
822 n3 = SIZE(basis%cm, 3)
823 ALLOCATE (gbasis%cm(n1, n2, 0:n3 - 1))
824 gbasis%cm = basis%cm
825 END IF
826 IF (ASSOCIATED(basis%as)) THEN
827 n1 = SIZE(basis%as, 1)
828 n2 = SIZE(basis%as, 2)
829 ALLOCATE (gbasis%as(n1, 0:n2 - 1))
830 gbasis%as = basis%as
831 END IF
832 IF (ASSOCIATED(basis%ns)) THEN
833 n1 = SIZE(basis%ns, 1)
834 n2 = SIZE(basis%ns, 2)
835 ALLOCATE (gbasis%ns(n1, 0:n2 - 1))
836 gbasis%ns = basis%ns
837 END IF
838 gbasis%eps_eig = basis%eps_eig
839 gbasis%geometrical = basis%geometrical
840 gbasis%aval = basis%aval
841 gbasis%cval = basis%cval
842 gbasis%start(0:lmat) = basis%start(0:lmat)
843
844 ! get information on quadrature type and number of grid points
845 ! allocate and initialize the atomic grid
846 NULLIFY (gbasis%grid)
847 CALL allocate_grid_atom(gbasis%grid)
848 ngp = SIZE(r)
849 IF (ngp <= 0) THEN
850 cpabort("The number of radial grid points must be greater than zero.")
851 END IF
852 cpassert(SIZE(rab) == ngp)
853 ! The supplied radial mesh needs no generated quadrature or angular grid.
854 ALLOCATE (gbasis%grid%rad(ngp), gbasis%grid%rad2(ngp), gbasis%grid%wr(ngp))
855 gbasis%grid%nr = ngp
856 gbasis%grid%rad(:) = r(:)
857 gbasis%grid%rad2(:) = r(:)*r(:)
858 gbasis%grid%wr(:) = rab(:)*gbasis%grid%rad2(:)
859
860 ! initialize basis function on a radial grid
861 nr = gbasis%grid%nr
862 m = maxval(gbasis%nbas)
863 ALLOCATE (gbasis%bf(nr, m, 0:lmat))
864 ALLOCATE (gbasis%dbf(nr, m, 0:lmat))
865 ALLOCATE (gbasis%ddbf(nr, m, 0:lmat))
866 gbasis%bf = 0._dp
867 gbasis%dbf = 0._dp
868 gbasis%ddbf = 0._dp
869
870 SELECT CASE (gbasis%basis_type)
871 CASE (gto_basis)
872 ! Zero prefactors must not evaluate negative powers at the UPF origin.
873 DO l = 0, lmat
874 DO i = 1, gbasis%nbas(l)
875 al = gbasis%am(i, l)
876 DO k = 1, nr
877 rk = gbasis%grid%rad(k)
878 ear = exp(-al*gbasis%grid%rad(k)**2)
879 gbasis%bf(k, i, l) = rk**l*ear
880 gbasis%dbf(k, i, l) = (real(l, dp)*rk**max(0, l - 1) - 2._dp*al*rk**(l + 1))*ear
881 gbasis%ddbf(k, i, l) = (real(l*(l - 1), dp)*rk**max(0, l - 2) - &
882 2._dp*al*real(2*l + 1, dp)*rk**(l) + 4._dp*al**2*rk**(l + 2))*ear
883 END DO
884 END DO
885 END DO
886 CASE (cgto_basis)
887 DO l = 0, lmat
888 DO i = 1, gbasis%nprim(l)
889 al = gbasis%am(i, l)
890 DO k = 1, nr
891 rk = gbasis%grid%rad(k)
892 ear = exp(-al*gbasis%grid%rad(k)**2)
893 DO j = 1, gbasis%nbas(l)
894 gbasis%bf(k, j, l) = gbasis%bf(k, j, l) + rk**l*ear*gbasis%cm(i, j, l)
895 gbasis%dbf(k, j, l) = gbasis%dbf(k, j, l) &
896 + (real(l, dp)*rk**max(0, l - 1) - 2._dp*al*rk**(l + 1))*ear*gbasis%cm(i, j, l)
897 gbasis%ddbf(k, j, l) = gbasis%ddbf(k, j, l) + &
898 (real(l*(l - 1), dp)*rk**max(0, l - 2) - 2._dp*al*real(2*l + 1, dp)*rk**(l) + 4._dp*al**2*rk**(l + 2))* &
899 ear*gbasis%cm(i, j, l)
900 END DO
901 END DO
902 END DO
903 END DO
904 CASE (sto_basis)
905 DO l = 0, lmat
906 DO i = 1, gbasis%nbas(l)
907 al = gbasis%as(i, l)
908 nl = gbasis%ns(i, l)
909 pf = (2._dp*al)**nl*sqrt(2._dp*al/fac(2*nl))
910 DO k = 1, nr
911 rk = gbasis%grid%rad(k)
912 ear = rk**(nl - 1)*exp(-al*rk)
913 gbasis%bf(k, i, l) = pf*ear
914 gbasis%dbf(k, i, l) = pf*(real(nl - 1, dp)/rk - al)*ear
915 gbasis%ddbf(k, i, l) = pf*(real((nl - 2)*(nl - 1), dp)/rk/rk &
916 - al*real(2*(nl - 1), dp)/rk + al*al)*ear
917 END DO
918 END DO
919 END DO
920 CASE (num_basis)
921 gbasis%basis_type = num_basis
922 cpabort("Numerical basis set type not yet implemented.")
923 CASE DEFAULT
924 cpabort("Unknown basis set type specified. Check the code!")
925 END SELECT
926
927 END SUBROUTINE atom_basis_gridrep
928
929! **************************************************************************************************
930!> \brief ...
931!> \param basis ...
932! **************************************************************************************************
933 SUBROUTINE release_atom_basis(basis)
934 TYPE(atom_basis_type), INTENT(INOUT) :: basis
935
936 IF (ASSOCIATED(basis%am)) THEN
937 DEALLOCATE (basis%am)
938 END IF
939 IF (ASSOCIATED(basis%cm)) THEN
940 DEALLOCATE (basis%cm)
941 END IF
942 IF (ASSOCIATED(basis%as)) THEN
943 DEALLOCATE (basis%as)
944 END IF
945 IF (ASSOCIATED(basis%ns)) THEN
946 DEALLOCATE (basis%ns)
947 END IF
948 IF (ASSOCIATED(basis%bf)) THEN
949 DEALLOCATE (basis%bf)
950 END IF
951 IF (ASSOCIATED(basis%dbf)) THEN
952 DEALLOCATE (basis%dbf)
953 END IF
954 IF (ASSOCIATED(basis%ddbf)) THEN
955 DEALLOCATE (basis%ddbf)
956 END IF
957
958 CALL deallocate_grid_atom(basis%grid)
959
960 END SUBROUTINE release_atom_basis
961! **************************************************************************************************
962
963! **************************************************************************************************
964!> \brief ...
965!> \param atom ...
966! **************************************************************************************************
967 SUBROUTINE create_atom_type(atom)
968 TYPE(atom_type), POINTER :: atom
969
970 cpassert(.NOT. ASSOCIATED(atom))
971
972 ALLOCATE (atom)
973
974 NULLIFY (atom%zmp_section)
975 NULLIFY (atom%xc_section)
976 NULLIFY (atom%fmat)
977 atom%do_zmp = .false.
978 atom%doread = .false.
979 atom%read_vxc = .false.
980 atom%dm = .false.
981 atom%hfx_pot%scale_coulomb = 0.0_dp
982 atom%hfx_pot%scale_longrange = 0.0_dp
983 atom%hfx_pot%omega = 0.0_dp
984
985 END SUBROUTINE create_atom_type
986
987! **************************************************************************************************
988!> \brief ...
989!> \param atom ...
990! **************************************************************************************************
991 SUBROUTINE release_atom_type(atom)
992 TYPE(atom_type), POINTER :: atom
993
994 cpassert(ASSOCIATED(atom))
995
996 NULLIFY (atom%basis)
997 NULLIFY (atom%integrals)
998 IF (ASSOCIATED(atom%state)) THEN
999 DEALLOCATE (atom%state)
1000 END IF
1001 IF (ASSOCIATED(atom%orbitals)) THEN
1002 CALL release_atom_orbs(atom%orbitals)
1003 END IF
1004
1005 IF (ASSOCIATED(atom%fmat)) CALL release_opmat(atom%fmat)
1006
1007 DEALLOCATE (atom)
1008
1009 END SUBROUTINE release_atom_type
1010
1011! ZMP adding input variables in subroutine do_zmp,doread,read_vxc,method_type
1012! **************************************************************************************************
1013!> \brief ...
1014!> \param atom ...
1015!> \param basis ...
1016!> \param state ...
1017!> \param integrals ...
1018!> \param orbitals ...
1019!> \param potential ...
1020!> \param zcore ...
1021!> \param pp_calc ...
1022!> \param do_zmp ...
1023!> \param doread ...
1024!> \param read_vxc ...
1025!> \param method_type ...
1026!> \param relativistic ...
1027!> \param coulomb_integral_type ...
1028!> \param exchange_integral_type ...
1029!> \param fmat ...
1030! **************************************************************************************************
1031 SUBROUTINE set_atom(atom, basis, state, integrals, orbitals, potential, zcore, pp_calc, do_zmp, doread, &
1032 read_vxc, method_type, relativistic, coulomb_integral_type, exchange_integral_type, fmat)
1033 TYPE(atom_type), POINTER :: atom
1034 TYPE(atom_basis_type), OPTIONAL, POINTER :: basis
1035 TYPE(atom_state), OPTIONAL, POINTER :: state
1036 TYPE(atom_integrals), OPTIONAL, POINTER :: integrals
1037 TYPE(atom_orbitals), OPTIONAL, POINTER :: orbitals
1038 TYPE(atom_potential_type), OPTIONAL, POINTER :: potential
1039 INTEGER, INTENT(IN), OPTIONAL :: zcore
1040 LOGICAL, INTENT(IN), OPTIONAL :: pp_calc, do_zmp, doread, read_vxc
1041 INTEGER, INTENT(IN), OPTIONAL :: method_type, relativistic, &
1042 coulomb_integral_type, &
1043 exchange_integral_type
1044 TYPE(opmat_type), OPTIONAL, POINTER :: fmat
1045
1046 cpassert(ASSOCIATED(atom))
1047
1048 IF (PRESENT(basis)) atom%basis => basis
1049 IF (PRESENT(state)) atom%state => state
1050 IF (PRESENT(integrals)) atom%integrals => integrals
1051 IF (PRESENT(orbitals)) atom%orbitals => orbitals
1052 IF (PRESENT(potential)) atom%potential => potential
1053 IF (PRESENT(zcore)) atom%zcore = zcore
1054 IF (PRESENT(pp_calc)) atom%pp_calc = pp_calc
1055! ZMP assigning variable values if present
1056 IF (PRESENT(do_zmp)) atom%do_zmp = do_zmp
1057 IF (PRESENT(doread)) atom%doread = doread
1058 IF (PRESENT(read_vxc)) atom%read_vxc = read_vxc
1059
1060 IF (PRESENT(method_type)) atom%method_type = method_type
1061 IF (PRESENT(relativistic)) atom%relativistic = relativistic
1062 IF (PRESENT(coulomb_integral_type)) atom%coulomb_integral_type = coulomb_integral_type
1063 IF (PRESENT(exchange_integral_type)) atom%exchange_integral_type = exchange_integral_type
1064
1065 IF (PRESENT(fmat)) THEN
1066 IF (ASSOCIATED(atom%fmat)) CALL release_opmat(atom%fmat)
1067 atom%fmat => fmat
1068 END IF
1069
1070 END SUBROUTINE set_atom
1071
1072! **************************************************************************************************
1073!> \brief ...
1074!> \param orbs ...
1075!> \param mbas ...
1076!> \param mo ...
1077! **************************************************************************************************
1078 SUBROUTINE create_atom_orbs(orbs, mbas, mo)
1079 TYPE(atom_orbitals), POINTER :: orbs
1080 INTEGER, INTENT(IN) :: mbas, mo
1081
1082 cpassert(.NOT. ASSOCIATED(orbs))
1083
1084 ALLOCATE (orbs)
1085
1086 ALLOCATE (orbs%wfn(mbas, mo, 0:lmat), orbs%wfna(mbas, mo, 0:lmat), orbs%wfnb(mbas, mo, 0:lmat))
1087 orbs%wfn = 0._dp
1088 orbs%wfna = 0._dp
1089 orbs%wfnb = 0._dp
1090
1091 ALLOCATE (orbs%pmat(mbas, mbas, 0:lmat), orbs%pmata(mbas, mbas, 0:lmat), orbs%pmatb(mbas, mbas, 0:lmat))
1092 orbs%pmat = 0._dp
1093 orbs%pmata = 0._dp
1094 orbs%pmatb = 0._dp
1095
1096 ALLOCATE (orbs%ener(mo, 0:lmat), orbs%enera(mo, 0:lmat), orbs%enerb(mo, 0:lmat))
1097 orbs%ener = 0._dp
1098 orbs%enera = 0._dp
1099 orbs%enerb = 0._dp
1100
1101 ALLOCATE (orbs%refene(mo, 0:lmat, 2), orbs%refchg(mo, 0:lmat, 2), orbs%refnod(mo, 0:lmat, 2))
1102 orbs%refene = 0._dp
1103 orbs%refchg = 0._dp
1104 orbs%refnod = 0._dp
1105 ALLOCATE (orbs%wrefene(mo, 0:lmat, 2), orbs%wrefchg(mo, 0:lmat, 2), orbs%wrefnod(mo, 0:lmat, 2))
1106 orbs%wrefene = 0._dp
1107 orbs%wrefchg = 0._dp
1108 orbs%wrefnod = 0._dp
1109 ALLOCATE (orbs%crefene(mo, 0:lmat, 2), orbs%crefchg(mo, 0:lmat, 2), orbs%crefnod(mo, 0:lmat, 2))
1110 orbs%crefene = 0._dp
1111 orbs%crefchg = 0._dp
1112 orbs%crefnod = 0._dp
1113 ALLOCATE (orbs%rcmax(mo, 0:lmat, 2))
1114 orbs%rcmax = 0._dp
1115 ALLOCATE (orbs%wpsir0(mo, 2), orbs%tpsir0(mo, 2))
1116 orbs%wpsir0 = 0._dp
1117 orbs%tpsir0 = 0._dp
1118 ALLOCATE (orbs%reftype(mo, 0:lmat, 2))
1119 orbs%reftype = "XX"
1120
1121 END SUBROUTINE create_atom_orbs
1122
1123! **************************************************************************************************
1124!> \brief ...
1125!> \param orbs ...
1126! **************************************************************************************************
1127 SUBROUTINE release_atom_orbs(orbs)
1128 TYPE(atom_orbitals), POINTER :: orbs
1129
1130 cpassert(ASSOCIATED(orbs))
1131
1132 IF (ASSOCIATED(orbs%wfn)) THEN
1133 DEALLOCATE (orbs%wfn, orbs%wfna, orbs%wfnb)
1134 END IF
1135 IF (ASSOCIATED(orbs%pmat)) THEN
1136 DEALLOCATE (orbs%pmat, orbs%pmata, orbs%pmatb)
1137 END IF
1138 IF (ASSOCIATED(orbs%ener)) THEN
1139 DEALLOCATE (orbs%ener, orbs%enera, orbs%enerb)
1140 END IF
1141 IF (ASSOCIATED(orbs%refene)) THEN
1142 DEALLOCATE (orbs%refene)
1143 END IF
1144 IF (ASSOCIATED(orbs%refchg)) THEN
1145 DEALLOCATE (orbs%refchg)
1146 END IF
1147 IF (ASSOCIATED(orbs%refnod)) THEN
1148 DEALLOCATE (orbs%refnod)
1149 END IF
1150 IF (ASSOCIATED(orbs%wrefene)) THEN
1151 DEALLOCATE (orbs%wrefene)
1152 END IF
1153 IF (ASSOCIATED(orbs%wrefchg)) THEN
1154 DEALLOCATE (orbs%wrefchg)
1155 END IF
1156 IF (ASSOCIATED(orbs%wrefnod)) THEN
1157 DEALLOCATE (orbs%wrefnod)
1158 END IF
1159 IF (ASSOCIATED(orbs%crefene)) THEN
1160 DEALLOCATE (orbs%crefene)
1161 END IF
1162 IF (ASSOCIATED(orbs%crefchg)) THEN
1163 DEALLOCATE (orbs%crefchg)
1164 END IF
1165 IF (ASSOCIATED(orbs%crefnod)) THEN
1166 DEALLOCATE (orbs%crefnod)
1167 END IF
1168 IF (ASSOCIATED(orbs%rcmax)) THEN
1169 DEALLOCATE (orbs%rcmax)
1170 END IF
1171 IF (ASSOCIATED(orbs%wpsir0)) THEN
1172 DEALLOCATE (orbs%wpsir0)
1173 END IF
1174 IF (ASSOCIATED(orbs%tpsir0)) THEN
1175 DEALLOCATE (orbs%tpsir0)
1176 END IF
1177 IF (ASSOCIATED(orbs%reftype)) THEN
1178 DEALLOCATE (orbs%reftype)
1179 END IF
1180
1181 DEALLOCATE (orbs)
1182
1183 END SUBROUTINE release_atom_orbs
1184
1185! **************************************************************************************************
1186!> \brief ...
1187!> \param hf_frac ...
1188!> \param do_hfx ...
1189!> \param atom ...
1190!> \param xc_section ...
1191!> \param extype ...
1192! **************************************************************************************************
1193 SUBROUTINE setup_hf_section(hf_frac, do_hfx, atom, xc_section, extype)
1194 REAL(kind=dp), INTENT(OUT) :: hf_frac
1195 LOGICAL, INTENT(OUT) :: do_hfx
1196 TYPE(atom_type), INTENT(IN), POINTER :: atom
1197 TYPE(section_vals_type), POINTER :: xc_section
1198 INTEGER, INTENT(IN) :: extype
1199
1200 INTEGER :: i, j, nr, nu, pot_type
1201 REAL(kind=dp) :: scale_coulomb, scale_longrange
1202 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: abscissa, weights
1203 TYPE(section_vals_type), POINTER :: hf_sub_section, hfx_sections
1204
1205 hf_frac = 0._dp
1206 IF (ASSOCIATED(atom%xc_section)) THEN
1207 xc_section => atom%xc_section
1208 hfx_sections => section_vals_get_subs_vals(xc_section, "HF")
1209 CALL section_vals_get(hfx_sections, explicit=do_hfx)
1210
1211 ! If nothing has been set explicitly, assume a Coulomb potential
1212 atom%hfx_pot%scale_longrange = 0.0_dp
1213 atom%hfx_pot%scale_coulomb = 1.0_dp
1214
1215 IF (do_hfx) THEN
1216 CALL section_vals_val_get(hfx_sections, "FRACTION", r_val=hf_frac)
1217
1218 ! Get potential info
1219 hf_sub_section => section_vals_get_subs_vals(hfx_sections, "INTERACTION_POTENTIAL", i_rep_section=1)
1220 CALL section_vals_val_get(hf_sub_section, "POTENTIAL_TYPE", i_val=pot_type)
1221 CALL section_vals_val_get(hf_sub_section, "OMEGA", r_val=atom%hfx_pot%omega)
1222 CALL section_vals_val_get(hf_sub_section, "SCALE_COULOMB", r_val=scale_coulomb)
1223 CALL section_vals_val_get(hf_sub_section, "SCALE_LONGRANGE", r_val=scale_longrange)
1224
1225 ! Setup atomic hfx potential
1226 SELECT CASE (pot_type)
1227 CASE DEFAULT
1228 cpwarn("Potential not implemented, use Coulomb instead!")
1230 atom%hfx_pot%scale_longrange = 0.0_dp
1231 atom%hfx_pot%scale_coulomb = scale_coulomb
1232 CASE (do_potential_long)
1233 atom%hfx_pot%scale_coulomb = 0.0_dp
1234 atom%hfx_pot%scale_longrange = scale_longrange
1235 CASE (do_potential_short)
1236 atom%hfx_pot%scale_coulomb = 1.0_dp
1237 atom%hfx_pot%scale_longrange = -1.0_dp
1238 CASE (do_potential_mix_cl)
1239 atom%hfx_pot%scale_coulomb = scale_coulomb
1240 atom%hfx_pot%scale_longrange = scale_longrange
1241 END SELECT
1242 END IF
1243
1244 ! Check whether extype is supported
1245 IF (atom%hfx_pot%scale_longrange /= 0.0_dp .AND. extype /= do_numeric .AND. extype /= do_semi_analytic) THEN
1246 cpabort("Only numerical and semi-analytic lrHF exchange available!")
1247 END IF
1248
1249 IF (atom%hfx_pot%scale_longrange /= 0.0_dp .AND. extype == do_numeric .AND. .NOT. ALLOCATED(atom%hfx_pot%kernel)) THEN
1250 CALL cite_reference(limpanuparb2011)
1251
1252 IF (atom%hfx_pot%do_gh) THEN
1253 ! Setup kernel for Ewald operator
1254 ! Because of the high computational costs of its calculation, we precalculate it here
1255 ! Use Gauss-Hermite grid instead of the external grid
1256 ALLOCATE (weights(atom%hfx_pot%nr_gh), abscissa(atom%hfx_pot%nr_gh))
1257 CALL get_gauss_hermite_weights(abscissa, weights, atom%hfx_pot%nr_gh)
1258
1259 nr = atom%basis%grid%nr
1260 ALLOCATE (atom%hfx_pot%kernel(nr, atom%hfx_pot%nr_gh, 0:atom%state%maxl_calc + atom%state%maxl_occ))
1261 atom%hfx_pot%kernel = 0.0_dp
1262 DO nu = 0, atom%state%maxl_calc + atom%state%maxl_occ
1263 DO i = 1, atom%hfx_pot%nr_gh
1264 DO j = 1, nr
1265 atom%hfx_pot%kernel(j, i, nu) = bessel0(2.0_dp*atom%hfx_pot%omega &
1266 *abscissa(i)*atom%basis%grid%rad(j), nu)*sqrt(weights(i))
1267 END DO
1268 END DO
1269 END DO
1270 ELSE
1271 ! Setup kernel for Ewald operator
1272 ! Because of the high computational costs of its calculation, we precalculate it here
1273 ! Choose it symmetric to further reduce the costs
1274 nr = atom%basis%grid%nr
1275 ALLOCATE (atom%hfx_pot%kernel(nr, nr, 0:atom%state%maxl_calc + atom%state%maxl_occ))
1276 atom%hfx_pot%kernel = 0.0_dp
1277 DO nu = 0, atom%state%maxl_calc + atom%state%maxl_occ
1278 DO i = 1, nr
1279 DO j = 1, i
1280 atom%hfx_pot%kernel(j, i, nu) = bessel0(2.0_dp*atom%hfx_pot%omega &
1281 *atom%basis%grid%rad(i)*atom%basis%grid%rad(j), nu)
1282 END DO
1283 END DO
1284 END DO
1285 END IF
1286 END IF
1287 ELSE
1288 NULLIFY (xc_section)
1289 do_hfx = .false.
1290 END IF
1291
1292 END SUBROUTINE setup_hf_section
1293
1294! **************************************************************************************************
1295!> \brief ...
1296!> \param abscissa ...
1297!> \param weights ...
1298!> \param nn ...
1299! **************************************************************************************************
1300 SUBROUTINE get_gauss_hermite_weights(abscissa, weights, nn)
1301 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: abscissa, weights
1302 INTEGER, INTENT(IN) :: nn
1303
1304 INTEGER :: counter, ii, info, liwork, lwork
1305 INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
1306 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: diag, subdiag, work
1307 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: eigenvec
1308
1309 ! Setup matrix for Golub-Welsch-algorithm to determine roots and weights of Gauss-Hermite quadrature
1310 ! If necessary, one can setup matrices differently for other quadratures
1311 ALLOCATE (work(1), iwork(1), diag(2*nn), subdiag(2*nn - 1), eigenvec(2*nn, 2*nn))
1312 lwork = -1
1313 liwork = -1
1314 diag = 0.0_dp
1315 DO ii = 1, 2*nn - 1
1316 subdiag(ii) = sqrt(real(ii, kind=dp)/2.0_dp)
1317 END DO
1318
1319 ! Get correct size for working matrices
1320 CALL dstevd('V', 2*nn, diag, subdiag, eigenvec, 2*nn, work, lwork, iwork, liwork, info)
1321 IF (info /= 0) THEN
1322 ! This should not happen!
1323 cpabort('Finding size of working matrices failed!')
1324 END IF
1325
1326 ! Setup working matrices with their respective optimal sizes
1327 lwork = int(work(1))
1328 liwork = iwork(1)
1329 DEALLOCATE (work, iwork)
1330 ALLOCATE (work(lwork), iwork(liwork))
1331
1332 ! Perform the actual eigenvalue decomposition
1333 CALL dstevd('V', 2*nn, diag, subdiag, eigenvec, 2*nn, work, lwork, iwork, liwork, info)
1334 IF (info /= 0) THEN
1335 ! This should not happen for the usual values of nn! (Checked for nn = 2000)
1336 cpabort('Eigenvalue decomposition failed!')
1337 END IF
1338
1339 DEALLOCATE (work, iwork, subdiag)
1340
1341 ! Identify positive roots of hermite polynomials (zeros of Hermite polynomials are symmetric wrt the origin)
1342 ! We will only keep the positive roots
1343 counter = 0
1344 DO ii = 1, 2*nn
1345 IF (diag(ii) > 0.0_dp) THEN
1346 counter = counter + 1
1347 abscissa(counter) = diag(ii)
1348 weights(counter) = rootpi*eigenvec(1, ii)**2
1349 END IF
1350 END DO
1351 IF (counter /= nn) THEN
1352 cpabort('Have not found enough or too many zeros!')
1353 END IF
1354
1355 END SUBROUTINE get_gauss_hermite_weights
1356
1357! **************************************************************************************************
1358!> \brief ...
1359!> \param opmat ...
1360!> \param n ...
1361!> \param lmax ...
1362! **************************************************************************************************
1363 SUBROUTINE create_opmat(opmat, n, lmax)
1364 TYPE(opmat_type), POINTER :: opmat
1365 INTEGER, DIMENSION(0:lmat), INTENT(IN) :: n
1366 INTEGER, INTENT(IN), OPTIONAL :: lmax
1367
1368 INTEGER :: lm, m
1369
1370 m = maxval(n)
1371 IF (PRESENT(lmax)) THEN
1372 lm = lmax
1373 ELSE
1374 lm = lmat
1375 END IF
1376
1377 cpassert(.NOT. ASSOCIATED(opmat))
1378
1379 ALLOCATE (opmat)
1380
1381 opmat%n = n
1382 ALLOCATE (opmat%op(m, m, 0:lm))
1383 opmat%op = 0._dp
1384
1385 END SUBROUTINE create_opmat
1386
1387! **************************************************************************************************
1388!> \brief ...
1389!> \param opmat ...
1390! **************************************************************************************************
1391 SUBROUTINE release_opmat(opmat)
1392 TYPE(opmat_type), POINTER :: opmat
1393
1394 cpassert(ASSOCIATED(opmat))
1395
1396 opmat%n = 0
1397 DEALLOCATE (opmat%op)
1398
1399 DEALLOCATE (opmat)
1400
1401 END SUBROUTINE release_opmat
1402
1403! **************************************************************************************************
1404!> \brief ...
1405!> \param opgrid ...
1406!> \param grid ...
1407! **************************************************************************************************
1408 SUBROUTINE create_opgrid(opgrid, grid)
1409 TYPE(opgrid_type), POINTER :: opgrid
1410 TYPE(grid_atom_type), POINTER :: grid
1411
1412 INTEGER :: nr
1413
1414 cpassert(.NOT. ASSOCIATED(opgrid))
1415
1416 ALLOCATE (opgrid)
1417
1418 opgrid%grid => grid
1419
1420 nr = grid%nr
1421
1422 ALLOCATE (opgrid%op(nr))
1423 opgrid%op = 0._dp
1424
1425 END SUBROUTINE create_opgrid
1426
1427! **************************************************************************************************
1428!> \brief ...
1429!> \param opgrid ...
1430! **************************************************************************************************
1431 SUBROUTINE release_opgrid(opgrid)
1432 TYPE(opgrid_type), POINTER :: opgrid
1433
1434 cpassert(ASSOCIATED(opgrid))
1435
1436 NULLIFY (opgrid%grid)
1437 DEALLOCATE (opgrid%op)
1438
1439 DEALLOCATE (opgrid)
1440
1441 END SUBROUTINE release_opgrid
1442
1443! **************************************************************************************************
1444!> \brief ...
1445!> \param zval ...
1446!> \param cval ...
1447!> \param aval ...
1448!> \param ngto ...
1449!> \param ival ...
1450! **************************************************************************************************
1451 SUBROUTINE clementi_geobas(zval, cval, aval, ngto, ival)
1452
1453 INTEGER, INTENT(IN) :: zval
1454 REAL(dp), INTENT(OUT) :: cval, aval
1455 INTEGER, DIMENSION(0:lmat), INTENT(OUT) :: ngto, ival
1456
1457 ngto = 0
1458 ival = 0
1459 cval = 0._dp
1460 aval = 0._dp
1461
1462 SELECT CASE (zval)
1463 CASE (1) ! this is from the general geometrical basis and extended
1464 cval = 2.0_dp
1465 aval = 0.016_dp
1466 ngto(0) = 20
1467 CASE (2)
1468 cval = 2.14774520_dp
1469 aval = 0.04850670_dp
1470 ngto(0) = 20
1471 CASE (3)
1472 cval = 2.08932430_dp
1473 aval = 0.02031060_dp
1474 ngto(0) = 23
1475 CASE (4)
1476 cval = 2.09753060_dp
1477 aval = 0.03207070_dp
1478 ngto(0) = 23
1479 CASE (5)
1480 cval = 2.10343410_dp
1481 aval = 0.03591970_dp
1482 ngto(0) = 23
1483 ngto(1) = 16
1484 CASE (6)
1485 cval = 2.10662820_dp
1486 aval = 0.05292410_dp
1487 ngto(0) = 23
1488 ngto(1) = 16
1489 CASE (7)
1490 cval = 2.13743840_dp
1491 aval = 0.06291970_dp
1492 ngto(0) = 23
1493 ngto(1) = 16
1494 CASE (8)
1495 cval = 2.08687310_dp
1496 aval = 0.08350860_dp
1497 ngto(0) = 23
1498 ngto(1) = 16
1499 CASE (9)
1500 cval = 2.12318180_dp
1501 aval = 0.09899170_dp
1502 ngto(0) = 23
1503 ngto(1) = 16
1504 CASE (10)
1505 cval = 2.13164810_dp
1506 aval = 0.11485350_dp
1507 ngto(0) = 23
1508 ngto(1) = 16
1509 CASE (11)
1510 cval = 2.11413310_dp
1511 aval = 0.00922630_dp
1512 ngto(0) = 26
1513 ngto(1) = 16
1514 ival(1) = 4
1515 CASE (12)
1516 cval = 2.12183620_dp
1517 aval = 0.01215850_dp
1518 ngto(0) = 26
1519 ngto(1) = 16
1520 ival(1) = 4
1521 CASE (13)
1522 cval = 2.06073230_dp
1523 aval = 0.01449350_dp
1524 ngto(0) = 26
1525 ngto(1) = 20
1526 ival(0) = 1
1527 CASE (14)
1528 cval = 2.08563660_dp
1529 aval = 0.01861460_dp
1530 ngto(0) = 26
1531 ngto(1) = 20
1532 ival(0) = 1
1533 CASE (15)
1534 cval = 2.04879270_dp
1535 aval = 0.02147790_dp
1536 ngto(0) = 26
1537 ngto(1) = 20
1538 ival(0) = 1
1539 CASE (16)
1540 cval = 2.06216660_dp
1541 aval = 0.01978920_dp
1542 ngto(0) = 26
1543 ngto(1) = 20
1544 ival(0) = 1
1545 CASE (17)
1546 cval = 2.04628670_dp
1547 aval = 0.02451470_dp
1548 ngto(0) = 26
1549 ngto(1) = 20
1550 ival(0) = 1
1551 CASE (18)
1552 cval = 2.08675200_dp
1553 aval = 0.02635040_dp
1554 ngto(0) = 26
1555 ngto(1) = 20
1556 ival(0) = 1
1557 CASE (19)
1558 cval = 2.02715220_dp
1559 aval = 0.01822040_dp
1560 ngto(0) = 29
1561 ngto(1) = 20
1562 ival(1) = 2
1563 CASE (20)
1564 cval = 2.01465650_dp
1565 aval = 0.01646570_dp
1566 ngto(0) = 29
1567 ngto(1) = 20
1568 ival(1) = 2
1569 CASE (21)
1570 cval = 2.01605240_dp
1571 aval = 0.01254190_dp
1572 ngto(0) = 30
1573 ngto(1) = 20
1574 ngto(2) = 18
1575 ival(1) = 2
1576 CASE (22)
1577 cval = 2.01800000_dp
1578 aval = 0.01195490_dp
1579 ngto(0) = 30
1580 ngto(1) = 21
1581 ngto(2) = 17
1582 ival(1) = 2
1583 ival(2) = 1
1584 CASE (23)
1585 cval = 1.98803560_dp
1586 aval = 0.02492140_dp
1587 ngto(0) = 30
1588 ngto(1) = 21
1589 ngto(2) = 17
1590 ival(1) = 2
1591 ival(2) = 1
1592 CASE (24)
1593 cval = 1.98984000_dp
1594 aval = 0.02568400_dp
1595 ngto(0) = 30
1596 ngto(1) = 21
1597 ngto(2) = 17
1598 ival(1) = 2
1599 ival(2) = 1
1600 CASE (25)
1601 cval = 2.01694380_dp
1602 aval = 0.02664480_dp
1603 ngto(0) = 30
1604 ngto(1) = 21
1605 ngto(2) = 17
1606 ival(1) = 2
1607 ival(2) = 1
1608 CASE (26)
1609 cval = 2.01824090_dp
1610 aval = 0.01355000_dp
1611 ngto(0) = 30
1612 ngto(1) = 21
1613 ngto(2) = 17
1614 ival(1) = 2
1615 ival(2) = 1
1616 CASE (27)
1617 cval = 1.98359400_dp
1618 aval = 0.01702210_dp
1619 ngto(0) = 30
1620 ngto(1) = 21
1621 ngto(2) = 17
1622 ival(1) = 2
1623 ival(2) = 2
1624 CASE (28)
1625 cval = 1.96797340_dp
1626 aval = 0.02163180_dp
1627 ngto(0) = 30
1628 ngto(1) = 22
1629 ngto(2) = 17
1630 ival(1) = 3
1631 ival(2) = 2
1632 CASE (29)
1633 cval = 1.98955180_dp
1634 aval = 0.02304480_dp
1635 ngto(0) = 30
1636 ngto(1) = 20
1637 ngto(2) = 17
1638 ival(1) = 3
1639 ival(2) = 2
1640 CASE (30)
1641 cval = 1.98074320_dp
1642 aval = 0.02754320_dp
1643 ngto(0) = 30
1644 ngto(1) = 21
1645 ngto(2) = 17
1646 ival(1) = 3
1647 ival(2) = 2
1648 CASE (31)
1649 cval = 2.00551070_dp
1650 aval = 0.02005530_dp
1651 ngto(0) = 30
1652 ngto(1) = 23
1653 ngto(2) = 17
1654 ival(0) = 1
1655 ival(2) = 2
1656 CASE (32)
1657 cval = 2.00000030_dp
1658 aval = 0.02003000_dp
1659 ngto(0) = 30
1660 ngto(1) = 24
1661 ngto(2) = 17
1662 ival(0) = 1
1663 ival(2) = 2
1664 CASE (33)
1665 cval = 2.00609100_dp
1666 aval = 0.02055620_dp
1667 ngto(0) = 30
1668 ngto(1) = 23
1669 ngto(2) = 17
1670 ival(0) = 1
1671 ival(2) = 2
1672 CASE (34)
1673 cval = 2.00701000_dp
1674 aval = 0.02230400_dp
1675 ngto(0) = 30
1676 ngto(1) = 24
1677 ngto(2) = 17
1678 ival(0) = 1
1679 ival(2) = 2
1680 CASE (35)
1681 cval = 2.01508710_dp
1682 aval = 0.02685790_dp
1683 ngto(0) = 30
1684 ngto(1) = 24
1685 ngto(2) = 17
1686 ival(0) = 1
1687 ival(2) = 2
1688 CASE (36)
1689 cval = 2.01960430_dp
1690 aval = 0.02960430_dp
1691 ngto(0) = 30
1692 ngto(1) = 24
1693 ngto(2) = 17
1694 ival(0) = 1
1695 ival(2) = 2
1696 CASE (37)
1697 cval = 2.00031000_dp
1698 aval = 0.00768400_dp
1699 ngto(0) = 32
1700 ngto(1) = 25
1701 ngto(2) = 17
1702 ival(0) = 1
1703 ival(1) = 1
1704 ival(2) = 4
1705 CASE (38)
1706 cval = 1.99563960_dp
1707 aval = 0.01401940_dp
1708 ngto(0) = 33
1709 ngto(1) = 24
1710 ngto(2) = 17
1711 ival(1) = 1
1712 ival(2) = 4
1713 CASE (39)
1714 cval = 1.98971210_dp
1715 aval = 0.01558470_dp
1716 ngto(0) = 33
1717 ngto(1) = 24
1718 ngto(2) = 20
1719 ival(1) = 1
1720 CASE (40)
1721 cval = 1.97976190_dp
1722 aval = 0.01705520_dp
1723 ngto(0) = 33
1724 ngto(1) = 24
1725 ngto(2) = 20
1726 ival(1) = 1
1727 CASE (41)
1728 cval = 1.97989290_dp
1729 aval = 0.01527040_dp
1730 ngto(0) = 33
1731 ngto(1) = 24
1732 ngto(2) = 20
1733 ival(1) = 1
1734 CASE (42)
1735 cval = 1.97909240_dp
1736 aval = 0.01879720_dp
1737 ngto(0) = 32
1738 ngto(1) = 24
1739 ngto(2) = 20
1740 ival(1) = 1
1741 CASE (43)
1742 cval = 1.98508430_dp
1743 aval = 0.01497550_dp
1744 ngto(0) = 32
1745 ngto(1) = 24
1746 ngto(2) = 20
1747 ival(1) = 2
1748 ival(2) = 1
1749 CASE (44)
1750 cval = 1.98515010_dp
1751 aval = 0.01856670_dp
1752 ngto(0) = 32
1753 ngto(1) = 24
1754 ngto(2) = 20
1755 ival(1) = 2
1756 ival(2) = 1
1757 CASE (45)
1758 cval = 1.98502970_dp
1759 aval = 0.01487000_dp
1760 ngto(0) = 32
1761 ngto(1) = 24
1762 ngto(2) = 20
1763 ival(1) = 2
1764 ival(2) = 1
1765 CASE (46)
1766 cval = 1.97672850_dp
1767 aval = 0.01762500_dp
1768 ngto(0) = 30
1769 ngto(1) = 24
1770 ngto(2) = 20
1771 ival(0) = 2
1772 ival(1) = 2
1773 ival(2) = 1
1774 CASE (47)
1775 cval = 1.97862730_dp
1776 aval = 0.01863310_dp
1777 ngto(0) = 32
1778 ngto(1) = 24
1779 ngto(2) = 20
1780 ival(1) = 2
1781 ival(2) = 1
1782 CASE (48)
1783 cval = 1.97990020_dp
1784 aval = 0.01347150_dp
1785 ngto(0) = 33
1786 ngto(1) = 24
1787 ngto(2) = 20
1788 ival(1) = 2
1789 ival(2) = 2
1790 CASE (49)
1791 cval = 1.97979410_dp
1792 aval = 0.00890265_dp
1793 ngto(0) = 33
1794 ngto(1) = 27
1795 ngto(2) = 20
1796 ival(0) = 2
1797 ival(2) = 2
1798 CASE (50)
1799 cval = 1.98001000_dp
1800 aval = 0.00895215_dp
1801 ngto(0) = 33
1802 ngto(1) = 27
1803 ngto(2) = 20
1804 ival(0) = 2
1805 ival(2) = 2
1806 CASE (51)
1807 cval = 1.97979980_dp
1808 aval = 0.01490290_dp
1809 ngto(0) = 33
1810 ngto(1) = 26
1811 ngto(2) = 20
1812 ival(1) = 1
1813 ival(2) = 2
1814 CASE (52)
1815 cval = 1.98009310_dp
1816 aval = 0.01490390_dp
1817 ngto(0) = 33
1818 ngto(1) = 26
1819 ngto(2) = 20
1820 ival(1) = 1
1821 ival(2) = 2
1822 CASE (53)
1823 cval = 1.97794750_dp
1824 aval = 0.01425880_dp
1825 ngto(0) = 33
1826 ngto(1) = 26
1827 ngto(2) = 20
1828 ival(0) = 2
1829 ival(1) = 1
1830 ival(2) = 2
1831 CASE (54)
1832 cval = 1.97784450_dp
1833 aval = 0.01430130_dp
1834 ngto(0) = 33
1835 ngto(1) = 26
1836 ngto(2) = 20
1837 ival(0) = 2
1838 ival(1) = 1
1839 ival(2) = 2
1840 CASE (55)
1841 cval = 1.97784450_dp
1842 aval = 0.00499318_dp
1843 ngto(0) = 32
1844 ngto(1) = 25
1845 ngto(2) = 17
1846 ival(0) = 1
1847 ival(1) = 3
1848 ival(2) = 6
1849 CASE (56)
1850 cval = 1.97764820_dp
1851 aval = 0.00500392_dp
1852 ngto(0) = 32
1853 ngto(1) = 25
1854 ngto(2) = 17
1855 ival(0) = 1
1856 ival(1) = 3
1857 ival(2) = 6
1858 CASE (57)
1859 cval = 1.97765150_dp
1860 aval = 0.00557083_dp
1861 ngto(0) = 32
1862 ngto(1) = 25
1863 ngto(2) = 20
1864 ival(0) = 1
1865 ival(1) = 3
1866 ival(2) = 3
1867 CASE (58)
1868 cval = 1.97768750_dp
1869 aval = 0.00547531_dp
1870 ngto(0) = 32
1871 ngto(1) = 25
1872 ngto(2) = 20
1873 ngto(3) = 16
1874 ival(0) = 1
1875 ival(1) = 3
1876 ival(2) = 3
1877 ival(3) = 3
1878 CASE (59)
1879 cval = 1.96986600_dp
1880 aval = 0.00813143_dp
1881 ngto(0) = 32
1882 ngto(1) = 25
1883 ngto(2) = 17
1884 ngto(3) = 16
1885 ival(0) = 1
1886 ival(1) = 3
1887 ival(2) = 6
1888 ival(3) = 4
1889 CASE (60)
1890 cval = 1.97765720_dp
1891 aval = 0.00489201_dp
1892 ngto(0) = 32
1893 ngto(1) = 25
1894 ngto(2) = 17
1895 ngto(3) = 16
1896 ival(0) = 1
1897 ival(1) = 3
1898 ival(2) = 6
1899 ival(3) = 4
1900 CASE (61)
1901 cval = 1.97768120_dp
1902 aval = 0.00499000_dp
1903 ngto(0) = 32
1904 ngto(1) = 25
1905 ngto(2) = 17
1906 ngto(3) = 16
1907 ival(0) = 1
1908 ival(1) = 3
1909 ival(2) = 6
1910 ival(3) = 4
1911 CASE (62)
1912 cval = 1.97745700_dp
1913 aval = 0.00615587_dp
1914 ngto(0) = 32
1915 ngto(1) = 25
1916 ngto(2) = 17
1917 ngto(3) = 16
1918 ival(0) = 1
1919 ival(1) = 3
1920 ival(2) = 6
1921 ival(3) = 4
1922 CASE (63)
1923 cval = 1.97570240_dp
1924 aval = 0.00769959_dp
1925 ngto(0) = 32
1926 ngto(1) = 25
1927 ngto(2) = 17
1928 ngto(3) = 16
1929 ival(0) = 1
1930 ival(1) = 3
1931 ival(2) = 6
1932 ival(3) = 4
1933 CASE (64)
1934 cval = 1.97629350_dp
1935 aval = 0.00706610_dp
1936 ngto(0) = 32
1937 ngto(1) = 25
1938 ngto(2) = 20
1939 ngto(3) = 16
1940 ival(0) = 1
1941 ival(1) = 3
1942 ival(2) = 3
1943 ival(3) = 4
1944 CASE (65)
1945 cval = 1.96900000_dp
1946 aval = 0.01019150_dp
1947 ngto(0) = 32
1948 ngto(1) = 26
1949 ngto(2) = 18
1950 ngto(3) = 16
1951 ival(0) = 1
1952 ival(1) = 3
1953 ival(2) = 6
1954 ival(3) = 4
1955 CASE (66)
1956 cval = 1.97350000_dp
1957 aval = 0.01334320_dp
1958 ngto(0) = 33
1959 ngto(1) = 26
1960 ngto(2) = 18
1961 ngto(3) = 16
1962 ival(0) = 1
1963 ival(1) = 3
1964 ival(2) = 6
1965 ival(3) = 4
1966 CASE (67)
1967 cval = 1.97493000_dp
1968 aval = 0.01331360_dp
1969 ngto(0) = 32
1970 ngto(1) = 24
1971 ngto(2) = 17
1972 ngto(3) = 14
1973 ival(1) = 2
1974 ival(2) = 5
1975 ival(3) = 4
1976 CASE (68)
1977 cval = 1.97597670_dp
1978 aval = 0.01434040_dp
1979 ngto(0) = 32
1980 ngto(1) = 24
1981 ngto(2) = 17
1982 ngto(3) = 14
1983 ival(0) = 0
1984 ival(1) = 2
1985 ival(2) = 5
1986 ival(3) = 4
1987 CASE (69)
1988 cval = 1.97809240_dp
1989 aval = 0.01529430_dp
1990 ngto(0) = 32
1991 ngto(1) = 24
1992 ngto(2) = 17
1993 ngto(3) = 14
1994 ival(0) = 0
1995 ival(1) = 2
1996 ival(2) = 5
1997 ival(3) = 4
1998 CASE (70)
1999 cval = 1.97644360_dp
2000 aval = 0.01312770_dp
2001 ngto(0) = 32
2002 ngto(1) = 24
2003 ngto(2) = 17
2004 ngto(3) = 14
2005 ival(0) = 0
2006 ival(1) = 2
2007 ival(2) = 5
2008 ival(3) = 4
2009 CASE (71)
2010 cval = 1.96998000_dp
2011 aval = 0.01745150_dp
2012 ngto(0) = 31
2013 ngto(1) = 24
2014 ngto(2) = 20
2015 ngto(3) = 14
2016 ival(0) = 1
2017 ival(1) = 2
2018 ival(2) = 2
2019 ival(3) = 4
2020 CASE (72)
2021 cval = 1.97223830_dp
2022 aval = 0.01639750_dp
2023 ngto(0) = 31
2024 ngto(1) = 24
2025 ngto(2) = 20
2026 ngto(3) = 14
2027 ival(0) = 1
2028 ival(1) = 2
2029 ival(2) = 2
2030 ival(3) = 4
2031 CASE (73)
2032 cval = 1.97462110_dp
2033 aval = 0.01603680_dp
2034 ngto(0) = 31
2035 ngto(1) = 24
2036 ngto(2) = 20
2037 ngto(3) = 14
2038 ival(0) = 1
2039 ival(1) = 2
2040 ival(2) = 2
2041 ival(3) = 4
2042 CASE (74)
2043 cval = 1.97756000_dp
2044 aval = 0.02030570_dp
2045 ngto(0) = 31
2046 ngto(1) = 24
2047 ngto(2) = 20
2048 ngto(3) = 14
2049 ival(0) = 1
2050 ival(1) = 2
2051 ival(2) = 2
2052 ival(3) = 4
2053 CASE (75)
2054 cval = 1.97645760_dp
2055 aval = 0.02057180_dp
2056 ngto(0) = 31
2057 ngto(1) = 24
2058 ngto(2) = 20
2059 ngto(3) = 14
2060 ival(0) = 1
2061 ival(1) = 2
2062 ival(2) = 2
2063 ival(3) = 4
2064 CASE (76)
2065 cval = 1.97725820_dp
2066 aval = 0.02058210_dp
2067 ngto(0) = 32
2068 ngto(1) = 24
2069 ngto(2) = 20
2070 ngto(3) = 15
2071 ival(0) = 0
2072 ival(1) = 2
2073 ival(2) = 2
2074 ival(3) = 4
2075 CASE (77)
2076 cval = 1.97749380_dp
2077 aval = 0.02219380_dp
2078 ngto(0) = 32
2079 ngto(1) = 24
2080 ngto(2) = 20
2081 ngto(3) = 15
2082 ival(0) = 0
2083 ival(1) = 2
2084 ival(2) = 2
2085 ival(3) = 4
2086 CASE (78)
2087 cval = 1.97946280_dp
2088 aval = 0.02216280_dp
2089 ngto(0) = 32
2090 ngto(1) = 24
2091 ngto(2) = 20
2092 ngto(3) = 15
2093 ival(0) = 0
2094 ival(1) = 2
2095 ival(2) = 2
2096 ival(3) = 4
2097 CASE (79)
2098 cval = 1.97852130_dp
2099 aval = 0.02168500_dp
2100 ngto(0) = 32
2101 ngto(1) = 24
2102 ngto(2) = 20
2103 ngto(3) = 15
2104 ival(0) = 0
2105 ival(1) = 2
2106 ival(2) = 2
2107 ival(3) = 4
2108 CASE (80)
2109 cval = 1.98045190_dp
2110 aval = 0.02177860_dp
2111 ngto(0) = 32
2112 ngto(1) = 24
2113 ngto(2) = 20
2114 ngto(3) = 15
2115 ival(0) = 0
2116 ival(1) = 2
2117 ival(2) = 2
2118 ival(3) = 4
2119 CASE (81)
2120 cval = 1.97000000_dp
2121 aval = 0.02275000_dp
2122 ngto(0) = 31
2123 ngto(1) = 25
2124 ngto(2) = 18
2125 ngto(3) = 13
2126 ival(0) = 1
2127 ival(1) = 0
2128 ival(2) = 3
2129 ival(3) = 6
2130 CASE (82)
2131 cval = 1.97713580_dp
2132 aval = 0.02317030_dp
2133 ngto(0) = 31
2134 ngto(1) = 27
2135 ngto(2) = 18
2136 ngto(3) = 13
2137 ival(0) = 1
2138 ival(1) = 0
2139 ival(2) = 3
2140 ival(3) = 6
2141 CASE (83)
2142 cval = 1.97537880_dp
2143 aval = 0.02672860_dp
2144 ngto(0) = 32
2145 ngto(1) = 27
2146 ngto(2) = 17
2147 ngto(3) = 13
2148 ival(0) = 1
2149 ival(1) = 0
2150 ival(2) = 3
2151 ival(3) = 6
2152 CASE (84)
2153 cval = 1.97545360_dp
2154 aval = 0.02745360_dp
2155 ngto(0) = 31
2156 ngto(1) = 27
2157 ngto(2) = 17
2158 ngto(3) = 13
2159 ival(0) = 1
2160 ival(1) = 0
2161 ival(2) = 3
2162 ival(3) = 6
2163 CASE (85)
2164 cval = 1.97338370_dp
2165 aval = 0.02616310_dp
2166 ngto(0) = 32
2167 ngto(1) = 27
2168 ngto(2) = 19
2169 ngto(3) = 13
2170 ival(0) = 1
2171 ival(1) = 0
2172 ival(2) = 3
2173 ival(3) = 6
2174 CASE (86)
2175 cval = 1.97294240_dp
2176 aval = 0.02429220_dp
2177 ngto(0) = 32
2178 ngto(1) = 27
2179 ngto(2) = 19
2180 ngto(3) = 13
2181 ival(0) = 1
2182 ival(1) = 0
2183 ival(2) = 3
2184 ival(3) = 6
2185 CASE (87:106) ! these numbers are an educated guess
2186 cval = 1.98000000_dp
2187 aval = 0.01400000_dp
2188 ngto(0) = 34
2189 ngto(1) = 28
2190 ngto(2) = 20
2191 ngto(3) = 15
2192 ival(0) = 0
2193 ival(1) = 0
2194 ival(2) = 3
2195 ival(3) = 6
2196 CASE DEFAULT
2197 cpabort("No geometrical basis set data are available for the selected atom number.")
2198 END SELECT
2199
2200 END SUBROUTINE clementi_geobas
2201
2202! **************************************************************************************************
2203!> \brief ...
2204!> \param element_symbol ...
2205!> \param basis ...
2206!> \param basis_set_name ...
2207!> \param basis_set_file ...
2208!> \param basis_section ...
2209! **************************************************************************************************
2210 SUBROUTINE read_basis_set(element_symbol, basis, basis_set_name, basis_set_file, &
2211 basis_section)
2212
2213 CHARACTER(LEN=*), INTENT(IN) :: element_symbol
2214 TYPE(atom_basis_type), INTENT(INOUT) :: basis
2215 CHARACTER(LEN=*), INTENT(IN) :: basis_set_name, basis_set_file
2216 TYPE(section_vals_type), POINTER :: basis_section
2217
2218 TYPE basis_set_type
2219 INTEGER :: lmin, lmax
2220 INTEGER, ALLOCATABLE :: l(:)
2221 REAL(dp), ALLOCATABLE :: zet(:), gcc(:, :)
2222 END TYPE basis_set_type
2223
2224 CHARACTER(len=20*default_string_length) :: line_att
2225 CHARACTER(LEN=240) :: line
2226 CHARACTER(LEN=242) :: line2
2227 CHARACTER(LEN=LEN(basis_set_name)) :: bsname
2228 CHARACTER(LEN=LEN(basis_set_name) + 2) :: bsname2
2229 CHARACTER(LEN=LEN(element_symbol)) :: symbol
2230 CHARACTER(LEN=LEN(element_symbol) + 2) :: symbol2
2231 INTEGER :: i, ii, ipgf, iset, ishell, j, k, &
2232 lshell, n, nj, npgf, ns, nset, nshell, &
2233 strlen1, strlen2
2234 INTEGER, ALLOCATABLE :: shells(:)
2235 LOGICAL :: found, is_ok, match, read_from_input
2236 REAL(dp) :: expzet, gcca, prefac, zeta
2237 TYPE(basis_set_type), ALLOCATABLE :: sets(:)
2238 TYPE(cp_parser_type) :: parser
2239 TYPE(cp_sll_val_type), POINTER :: list
2240 TYPE(val_type), POINTER :: val
2241
2242 bsname = basis_set_name
2243 symbol = element_symbol
2244
2245 read_from_input = .false.
2246 CALL section_vals_get(basis_section, explicit=read_from_input)
2247 IF (read_from_input) THEN
2248 NULLIFY (list, val)
2249 CALL section_vals_list_get(basis_section, "_DEFAULT_KEYWORD_", list=list)
2250 CALL uppercase(symbol)
2251 CALL uppercase(bsname)
2252 is_ok = cp_sll_val_next(list, val)
2253 cpassert(is_ok)
2254 CALL val_get(val, c_val=line_att)
2255 READ (line_att, *) nset
2256 ELSE
2257 CALL parser_create(parser, basis_set_file)
2258 ! Search for the requested basis set in the basis set file
2259 ! until the basis set is found or the end of file is reached
2260 search_loop: DO
2261 CALL parser_search_string(parser, trim(bsname), .true., found, line)
2262 IF (found) THEN
2263 CALL uppercase(symbol)
2264 CALL uppercase(bsname)
2265 match = .false.
2266 CALL uppercase(line)
2267 ! Check both the element symbol and the basis set name
2268 line2 = " "//line//" "
2269 symbol2 = " "//trim(symbol)//" "
2270 bsname2 = " "//trim(bsname)//" "
2271 strlen1 = len_trim(symbol2) + 1
2272 strlen2 = len_trim(bsname2) + 1
2273
2274 IF ((index(line2, symbol2(:strlen1)) > 0) .AND. &
2275 (index(line2, bsname2(:strlen2)) > 0)) match = .true.
2276
2277 IF (match) THEN
2278 ! Read the basis set information
2279 CALL parser_get_object(parser, nset, newline=.true.)
2280 EXIT search_loop
2281
2282 END IF
2283 ELSE
2284 ! Stop program, if the end of file is reached
2285 cpabort("End of file reached and the requested basis set was not found.")
2286 END IF
2287
2288 END DO search_loop
2289 END IF
2290
2291 cpassert(nset > 0)
2292 ALLOCATE (sets(nset))
2293 DO iset = 1, nset
2294 associate(set => sets(iset))
2295 IF (read_from_input) THEN
2296 is_ok = cp_sll_val_next(list, val)
2297 cpassert(is_ok)
2298 CALL val_get(val, c_val=line_att)
2299 READ (line_att, *) n, set%lmin, set%lmax, npgf
2300 DO i = 1, 4
2301 CALL remove_word(line_att)
2302 END DO
2303 ELSE
2304 CALL parser_get_object(parser, n, newline=.true.)
2305 CALL parser_get_object(parser, set%lmin)
2306 CALL parser_get_object(parser, set%lmax)
2307 CALL parser_get_object(parser, npgf)
2308 END IF
2309 cpassert(set%lmin >= 0 .AND. set%lmax >= set%lmin .AND. npgf > 0)
2310 ALLOCATE (shells(set%lmin:set%lmax))
2311 DO lshell = set%lmin, set%lmax
2312 IF (read_from_input) THEN
2313 READ (line_att, *) shells(lshell)
2314 CALL remove_word(line_att)
2315 ELSE
2316 CALL parser_get_object(parser, shells(lshell))
2317 END IF
2318 END DO
2319 IF (read_from_input) THEN
2320 cpassert(len_trim(line_att) == 0)
2321 END IF
2322 cpassert(all(shells >= 0))
2323 nshell = sum(shells)
2324 ALLOCATE (set%l(nshell), set%zet(npgf), set%gcc(npgf, nshell))
2325 ishell = 0
2326 DO lshell = set%lmin, set%lmax
2327 set%l(ishell + 1:ishell + shells(lshell)) = lshell
2328 ishell = ishell + shells(lshell)
2329 END DO
2330 DEALLOCATE (shells)
2331 DO ipgf = 1, npgf
2332 IF (read_from_input) THEN
2333 is_ok = cp_sll_val_next(list, val)
2334 cpassert(is_ok)
2335 CALL val_get(val, c_val=line_att)
2336 READ (line_att, *) set%zet(ipgf), set%gcc(ipgf, :)
2337 ELSE
2338 CALL parser_get_object(parser, set%zet(ipgf), newline=.true.)
2339 DO ishell = 1, nshell
2340 CALL parser_get_object(parser, set%gcc(ipgf, ishell))
2341 END DO
2342 END IF
2343 END DO
2344 END associate
2345 END DO
2346 IF (.NOT. read_from_input) CALL parser_release(parser)
2347
2348 ! fill in the basis data structures
2349 basis%nprim = 0
2350 basis%nbas = 0
2351 DO i = 1, nset
2352 DO j = sets(i)%lmin, min(sets(i)%lmax, lmat)
2353 basis%nprim(j) = basis%nprim(j) + SIZE(sets(i)%zet)
2354 END DO
2355 DO j = 1, SIZE(sets(i)%l)
2356 k = sets(i)%l(j)
2357 IF (k <= lmat) basis%nbas(k) = basis%nbas(k) + 1
2358 END DO
2359 END DO
2360
2361 nj = maxval(basis%nprim)
2362 ns = maxval(basis%nbas)
2363 ALLOCATE (basis%am(nj, 0:lmat))
2364 basis%am = 0._dp
2365 ALLOCATE (basis%cm(nj, ns, 0:lmat))
2366 basis%cm = 0._dp
2367
2368 DO j = 0, lmat
2369 nj = 0
2370 ns = 0
2371 DO i = 1, nset
2372 IF (j >= sets(i)%lmin .AND. j <= sets(i)%lmax) THEN
2373 DO ipgf = 1, SIZE(sets(i)%zet)
2374 basis%am(nj + ipgf, j) = sets(i)%zet(ipgf)
2375 END DO
2376 DO ii = 1, SIZE(sets(i)%l)
2377 IF (sets(i)%l(ii) == j) THEN
2378 ns = ns + 1
2379 DO ipgf = 1, SIZE(sets(i)%zet)
2380 basis%cm(nj + ipgf, ns, j) = sets(i)%gcc(ipgf, ii)
2381 END DO
2382 END IF
2383 END DO
2384 nj = nj + SIZE(sets(i)%zet)
2385 END IF
2386 END DO
2387 END DO
2388
2389 ! Normalization
2390 DO j = 0, lmat
2391 expzet = 0.25_dp*real(2*j + 3, dp)
2392 prefac = sqrt(rootpi/2._dp**(j + 2)*dfac(2*j + 1))
2393 DO ipgf = 1, basis%nprim(j)
2394 DO ii = 1, basis%nbas(j)
2395 gcca = basis%cm(ipgf, ii, j)
2396 zeta = 2._dp*basis%am(ipgf, j)
2397 basis%cm(ipgf, ii, j) = zeta**expzet*gcca/prefac
2398 END DO
2399 END DO
2400 END DO
2401
2402 END SUBROUTINE read_basis_set
2403
2404! **************************************************************************************************
2405!> \brief ...
2406!> \param optimization ...
2407!> \param opt_section ...
2408! **************************************************************************************************
2409 SUBROUTINE read_atom_opt_section(optimization, opt_section)
2410 TYPE(atom_optimization_type), INTENT(INOUT) :: optimization
2411 TYPE(section_vals_type), POINTER :: opt_section
2412
2413 INTEGER :: miter, ndiis
2414 REAL(kind=dp) :: damp, eps_diis, eps_scf
2415
2416 CALL section_vals_val_get(opt_section, "MAX_ITER", i_val=miter)
2417 CALL section_vals_val_get(opt_section, "EPS_SCF", r_val=eps_scf)
2418 CALL section_vals_val_get(opt_section, "N_DIIS", i_val=ndiis)
2419 CALL section_vals_val_get(opt_section, "EPS_DIIS", r_val=eps_diis)
2420 CALL section_vals_val_get(opt_section, "DAMPING", r_val=damp)
2421
2422 optimization%max_iter = miter
2423 optimization%eps_scf = eps_scf
2424 optimization%n_diis = ndiis
2425 optimization%eps_diis = eps_diis
2426 optimization%damping = damp
2427
2428 END SUBROUTINE read_atom_opt_section
2429! **************************************************************************************************
2430!> \brief ...
2431!> \param potential ...
2432!> \param potential_section ...
2433!> \param zval ...
2434! **************************************************************************************************
2435 SUBROUTINE init_atom_potential(potential, potential_section, zval)
2436 TYPE(atom_potential_type), INTENT(INOUT) :: potential
2437 TYPE(section_vals_type), POINTER :: potential_section
2438 INTEGER, INTENT(IN) :: zval
2439
2440 CHARACTER(LEN=default_path_length) :: pseudo_fn
2441 CHARACTER(LEN=default_string_length) :: pseudo_name
2442 INTEGER :: ic
2443 REAL(dp), DIMENSION(:), POINTER :: convals
2444 TYPE(section_vals_type), POINTER :: ecp_potential_section, &
2445 gth_potential_section
2446
2447 IF (zval > 0) THEN
2448 CALL section_vals_val_get(potential_section, "PSEUDO_TYPE", i_val=potential%ppot_type)
2449
2450 SELECT CASE (potential%ppot_type)
2451 CASE (gth_pseudo)
2452 CALL section_vals_val_get(potential_section, "POTENTIAL_FILE_NAME", c_val=pseudo_fn)
2453 CALL section_vals_val_get(potential_section, "POTENTIAL_NAME", c_val=pseudo_name)
2454 gth_potential_section => section_vals_get_subs_vals(potential_section, "GTH_POTENTIAL")
2455 CALL read_gth_potential(ptable(zval)%symbol, potential%gth_pot, &
2456 pseudo_name, pseudo_fn, gth_potential_section)
2457 CASE (ecp_pseudo)
2458 CALL section_vals_val_get(potential_section, "POTENTIAL_FILE_NAME", c_val=pseudo_fn)
2459 CALL section_vals_val_get(potential_section, "POTENTIAL_NAME", c_val=pseudo_name)
2460 ecp_potential_section => section_vals_get_subs_vals(potential_section, "ECP")
2461 CALL read_ecp_potential(ptable(zval)%symbol, potential%ecp_pot, &
2462 pseudo_name, pseudo_fn, ecp_potential_section)
2463 CASE (upf_pseudo)
2464 CALL section_vals_val_get(potential_section, "POTENTIAL_FILE_NAME", c_val=pseudo_fn)
2465 CALL section_vals_val_get(potential_section, "POTENTIAL_NAME", c_val=pseudo_name)
2466 CALL atom_read_upf(potential%upf_pot, pseudo_fn)
2467 potential%upf_pot%pname = pseudo_name
2468 CASE (sgp_pseudo)
2469 cpabort("Pseudopotential type SGP is not implemented.")
2470 CASE (no_pseudo)
2471 ! do nothing
2472 CASE DEFAULT
2473 cpabort("Invalid pseudopotential type selected. Check the code!")
2474 END SELECT
2475 ELSE
2476 potential%ppot_type = no_pseudo
2477 END IF
2478
2479 ! confinement
2480 NULLIFY (convals)
2481 CALL section_vals_val_get(potential_section, "CONFINEMENT_TYPE", i_val=ic)
2482 potential%conf_type = ic
2483 IF (potential%conf_type == no_conf) THEN
2484 potential%acon = 0.0_dp
2485 potential%rcon = 4.0_dp
2486 potential%scon = 2.0_dp
2487 potential%confinement = .false.
2488 ELSE IF (potential%conf_type == poly_conf) THEN
2489 CALL section_vals_val_get(potential_section, "CONFINEMENT", r_vals=convals)
2490 IF (SIZE(convals) >= 1) THEN
2491 IF (convals(1) > 0.0_dp) THEN
2492 potential%confinement = .true.
2493 potential%acon = convals(1)
2494 IF (SIZE(convals) >= 2) THEN
2495 potential%rcon = convals(2)
2496 ELSE
2497 potential%rcon = 4.0_dp
2498 END IF
2499 IF (SIZE(convals) >= 3) THEN
2500 potential%scon = convals(3)
2501 ELSE
2502 potential%scon = 2.0_dp
2503 END IF
2504 ELSE
2505 potential%confinement = .false.
2506 END IF
2507 ELSE
2508 potential%confinement = .false.
2509 END IF
2510 ELSE IF (potential%conf_type == barrier_conf) THEN
2511 potential%acon = 200.0_dp
2512 potential%rcon = 4.0_dp
2513 potential%scon = 12.0_dp
2514 potential%confinement = .true.
2515 CALL section_vals_val_get(potential_section, "CONFINEMENT", r_vals=convals)
2516 IF (SIZE(convals) >= 1) THEN
2517 IF (convals(1) > 0.0_dp) THEN
2518 potential%acon = convals(1)
2519 IF (SIZE(convals) >= 2) THEN
2520 potential%rcon = convals(2)
2521 END IF
2522 IF (SIZE(convals) >= 3) THEN
2523 potential%scon = convals(3)
2524 END IF
2525 ELSE
2526 potential%confinement = .false.
2527 END IF
2528 END IF
2529 END IF
2530
2531 END SUBROUTINE init_atom_potential
2532! **************************************************************************************************
2533!> \brief ...
2534!> \param potential ...
2535! **************************************************************************************************
2536 SUBROUTINE release_atom_potential(potential)
2537 TYPE(atom_potential_type), INTENT(INOUT) :: potential
2538
2539 potential%confinement = .false.
2540
2541 CALL atom_release_upf(potential%upf_pot)
2542
2543 END SUBROUTINE release_atom_potential
2544! **************************************************************************************************
2545!> \brief ...
2546!> \param element_symbol ...
2547!> \param potential ...
2548!> \param pseudo_name ...
2549!> \param pseudo_file ...
2550!> \param potential_section ...
2551! **************************************************************************************************
2552 SUBROUTINE read_gth_potential(element_symbol, potential, pseudo_name, pseudo_file, &
2553 potential_section)
2554
2555 CHARACTER(LEN=*), INTENT(IN) :: element_symbol
2556 TYPE(atom_gthpot_type), INTENT(INOUT) :: potential
2557 CHARACTER(LEN=*), INTENT(IN) :: pseudo_name, pseudo_file
2558 TYPE(section_vals_type), POINTER :: potential_section
2559
2560 CHARACTER(LEN=240) :: line
2561 CHARACTER(LEN=242) :: line2
2562 CHARACTER(len=5*default_string_length) :: line_att
2563 CHARACTER(LEN=LEN(element_symbol)) :: symbol
2564 CHARACTER(LEN=LEN(element_symbol)+2) :: symbol2
2565 CHARACTER(LEN=LEN(pseudo_name)) :: apname
2566 CHARACTER(LEN=LEN(pseudo_name)+2) :: apname2
2567 INTEGER :: i, ic, ipot, j, l, nlmax, strlen1, &
2568 strlen2
2569 INTEGER, DIMENSION(0:lmat) :: elec_conf
2570 LOGICAL :: found, is_ok, match, read_from_input
2571 TYPE(cp_sll_val_type), POINTER :: list
2572 TYPE(val_type), POINTER :: val
2573
2574 elec_conf = 0
2575
2576 apname = pseudo_name
2577 symbol = element_symbol
2578
2579 potential%symbol = symbol
2580 potential%pname = apname
2581 potential%econf = 0
2582 potential%rc = 0._dp
2583 potential%ncl = 0
2584 potential%cl = 0._dp
2585 potential%nl = 0
2586 potential%rcnl = 0._dp
2587 potential%hnl = 0._dp
2588 potential%soc = .false.
2589 potential%knl = 0._dp
2590
2591 potential%lpotextended = .false.
2592 potential%lsdpot = .false.
2593 potential%nlcc = .false.
2594 potential%nexp_lpot = 0
2595 potential%nexp_lsd = 0
2596 potential%nexp_nlcc = 0
2597
2598 read_from_input = .false.
2599 CALL section_vals_get(potential_section, explicit=read_from_input)
2600 IF (read_from_input) THEN
2601 CALL section_vals_list_get(potential_section, "_DEFAULT_KEYWORD_", list=list)
2602 CALL uppercase(symbol)
2603 CALL uppercase(apname)
2604 ! Read the electronic configuration, not used here
2605 l = 0
2606 is_ok = cp_sll_val_next(list, val)
2607 cpassert(is_ok)
2608 CALL val_get(val, c_val=line_att)
2609 READ (line_att, *) elec_conf(l)
2610 CALL remove_word(line_att)
2611 DO WHILE (len_trim(line_att) /= 0)
2612 l = l + 1
2613 READ (line_att, *) elec_conf(l)
2614 CALL remove_word(line_att)
2615 END DO
2616 potential%econf(0:lmat) = elec_conf(0:lmat)
2617 potential%zion = real(sum(elec_conf), dp)
2618 ! Read r(loc) to define the exponent of the core charge
2619 is_ok = cp_sll_val_next(list, val)
2620 cpassert(is_ok)
2621 CALL val_get(val, c_val=line_att)
2622 READ (line_att, *) potential%rc
2623 CALL remove_word(line_att)
2624 ! Read the parameters for the local part of the GTH pseudopotential (ppl)
2625 READ (line_att, *) potential%ncl
2626 CALL remove_word(line_att)
2627 DO i = 1, potential%ncl
2628 READ (line_att, *) potential%cl(i)
2629 CALL remove_word(line_att)
2630 END DO
2631 ! Check for the next entry: LPOT, NLCC, LSD, or ppnl
2632 DO
2633 is_ok = cp_sll_val_next(list, val)
2634 cpassert(is_ok)
2635 CALL val_get(val, c_val=line_att)
2636 IF (index(line_att, "LPOT") /= 0) THEN
2637 potential%lpotextended = .true.
2638 CALL remove_word(line_att)
2639 READ (line_att, *) potential%nexp_lpot
2640 DO ipot = 1, potential%nexp_lpot
2641 is_ok = cp_sll_val_next(list, val)
2642 cpassert(is_ok)
2643 CALL val_get(val, c_val=line_att)
2644 READ (line_att, *) potential%alpha_lpot(ipot)
2645 CALL remove_word(line_att)
2646 READ (line_att, *) potential%nct_lpot(ipot)
2647 CALL remove_word(line_att)
2648 DO ic = 1, potential%nct_lpot(ipot)
2649 READ (line_att, *) potential%cval_lpot(ic, ipot)
2650 CALL remove_word(line_att)
2651 END DO
2652 END DO
2653 ELSE IF (index(line_att, "NLCC") /= 0) THEN
2654 potential%nlcc = .true.
2655 CALL remove_word(line_att)
2656 READ (line_att, *) potential%nexp_nlcc
2657 DO ipot = 1, potential%nexp_nlcc
2658 is_ok = cp_sll_val_next(list, val)
2659 cpassert(is_ok)
2660 CALL val_get(val, c_val=line_att)
2661 READ (line_att, *) potential%alpha_nlcc(ipot)
2662 CALL remove_word(line_att)
2663 READ (line_att, *) potential%nct_nlcc(ipot)
2664 CALL remove_word(line_att)
2665 DO ic = 1, potential%nct_nlcc(ipot)
2666 READ (line_att, *) potential%cval_nlcc(ic, ipot)
2667 !make cp2k compatible with bigdft
2668 potential%cval_nlcc(ic, ipot) = potential%cval_nlcc(ic, ipot)/(4.0_dp*pi)
2669 CALL remove_word(line_att)
2670 END DO
2671 END DO
2672 ELSE IF (index(line_att, "LSD") /= 0) THEN
2673 potential%lsdpot = .true.
2674 CALL remove_word(line_att)
2675 READ (line_att, *) potential%nexp_lsd
2676 DO ipot = 1, potential%nexp_lsd
2677 is_ok = cp_sll_val_next(list, val)
2678 cpassert(is_ok)
2679 CALL val_get(val, c_val=line_att)
2680 READ (line_att, *) potential%alpha_lsd(ipot)
2681 CALL remove_word(line_att)
2682 READ (line_att, *) potential%nct_lsd(ipot)
2683 CALL remove_word(line_att)
2684 DO ic = 1, potential%nct_lsd(ipot)
2685 READ (line_att, *) potential%cval_lsd(ic, ipot)
2686 CALL remove_word(line_att)
2687 END DO
2688 END DO
2689 ELSE
2690 EXIT
2691 END IF
2692 END DO
2693 ! Read the parameters for the non-local part of the GTH pseudopotential (ppnl)
2694 READ (line_att, *) nlmax
2695 CALL remove_word(line_att)
2696 IF (index(line_att, "SOC") /= 0) potential%soc = .true.
2697 IF (nlmax > 0) THEN
2698 ! Load the parameter for nlmax non-local projectors
2699 DO l = 0, nlmax - 1
2700 is_ok = cp_sll_val_next(list, val)
2701 cpassert(is_ok)
2702 CALL val_get(val, c_val=line_att)
2703 READ (line_att, *) potential%rcnl(l)
2704 CALL remove_word(line_att)
2705 READ (line_att, *) potential%nl(l)
2706 CALL remove_word(line_att)
2707 DO i = 1, potential%nl(l)
2708 IF (i == 1) THEN
2709 READ (line_att, *) potential%hnl(1, 1, l)
2710 CALL remove_word(line_att)
2711 ELSE
2712 cpassert(len_trim(line_att) == 0)
2713 is_ok = cp_sll_val_next(list, val)
2714 cpassert(is_ok)
2715 CALL val_get(val, c_val=line_att)
2716 READ (line_att, *) potential%hnl(i, i, l)
2717 CALL remove_word(line_att)
2718 END IF
2719 DO j = i + 1, potential%nl(l)
2720 READ (line_att, *) potential%hnl(i, j, l)
2721 potential%hnl(j, i, l) = potential%hnl(i, j, l)
2722 CALL remove_word(line_att)
2723 END DO
2724 END DO
2725 IF (potential%soc .AND. l /= 0) THEN
2726 is_ok = cp_sll_val_next(list, val)
2727 cpassert(is_ok)
2728 CALL val_get(val, c_val=line_att)
2729 DO i = 1, potential%nl(l)
2730 IF (i == 1) THEN
2731 READ (line_att, *) potential%knl(1, 1, l)
2732 CALL remove_word(line_att)
2733 ELSE
2734 cpassert(len_trim(line_att) == 0)
2735 is_ok = cp_sll_val_next(list, val)
2736 cpassert(is_ok)
2737 CALL val_get(val, c_val=line_att)
2738 READ (line_att, *) potential%knl(i, i, l)
2739 CALL remove_word(line_att)
2740 END IF
2741 DO j = i + 1, potential%nl(l)
2742 READ (line_att, *) potential%knl(i, j, l)
2743 potential%knl(j, i, l) = potential%knl(i, j, l)
2744 CALL remove_word(line_att)
2745 END DO
2746 END DO
2747 END IF
2748 cpassert(len_trim(line_att) == 0)
2749 END DO
2750 END IF
2751 ELSE
2752 block
2753 TYPE(cp_parser_type) :: parser
2754 CALL parser_create(parser, pseudo_file)
2755
2756 search_loop: DO
2757 CALL parser_search_string(parser, trim(apname), .true., found, line)
2758 IF (found) THEN
2759 CALL uppercase(symbol)
2760 CALL uppercase(apname)
2761 ! Check both the element symbol and the atomic potential name
2762 match = .false.
2763 CALL uppercase(line)
2764 line2 = " "//line//" "
2765 symbol2 = " "//trim(symbol)//" "
2766 apname2 = " "//trim(apname)//" "
2767 strlen1 = len_trim(symbol2) + 1
2768 strlen2 = len_trim(apname2) + 1
2769
2770 IF ((index(line2, symbol2(:strlen1)) > 0) .AND. &
2771 (index(line2, apname2(:strlen2)) > 0)) match = .true.
2772
2773 IF (match) THEN
2774 ! Read the electronic configuration
2775 l = 0
2776 CALL parser_get_object(parser, elec_conf(l), newline=.true.)
2777 DO WHILE (parser_test_next_token(parser) == "INT")
2778 l = l + 1
2779 CALL parser_get_object(parser, elec_conf(l))
2780 END DO
2781 potential%econf(0:lmat) = elec_conf(0:lmat)
2782 potential%zion = real(sum(elec_conf), dp)
2783 ! Read r(loc) to define the exponent of the core charge
2784 CALL parser_get_object(parser, potential%rc, newline=.true.)
2785 ! Read the parameters for the local part of the GTH pseudopotential (ppl)
2786 CALL parser_get_object(parser, potential%ncl)
2787 DO i = 1, potential%ncl
2788 CALL parser_get_object(parser, potential%cl(i))
2789 END DO
2790 ! Extended type input
2791 DO
2792 CALL parser_get_next_line(parser, 1)
2793 IF (parser_test_next_token(parser) == "INT") THEN
2794 EXIT
2795 ELSE IF (parser_test_next_token(parser) == "STR") THEN
2796 CALL parser_get_object(parser, line)
2797 IF (index(line, "LPOT") /= 0) THEN
2798 ! local potential
2799 potential%lpotextended = .true.
2800 CALL parser_get_object(parser, potential%nexp_lpot)
2801 DO ipot = 1, potential%nexp_lpot
2802 CALL parser_get_object(parser, potential%alpha_lpot(ipot), newline=.true.)
2803 CALL parser_get_object(parser, potential%nct_lpot(ipot))
2804 DO ic = 1, potential%nct_lpot(ipot)
2805 CALL parser_get_object(parser, potential%cval_lpot(ic, ipot))
2806 END DO
2807 END DO
2808 ELSE IF (index(line, "NLCC") /= 0) THEN
2809 ! NLCC
2810 potential%nlcc = .true.
2811 CALL parser_get_object(parser, potential%nexp_nlcc)
2812 DO ipot = 1, potential%nexp_nlcc
2813 CALL parser_get_object(parser, potential%alpha_nlcc(ipot), newline=.true.)
2814 CALL parser_get_object(parser, potential%nct_nlcc(ipot))
2815 DO ic = 1, potential%nct_nlcc(ipot)
2816 CALL parser_get_object(parser, potential%cval_nlcc(ic, ipot))
2817 !make cp2k compatible with bigdft
2818 potential%cval_nlcc(ic, ipot) = potential%cval_nlcc(ic, ipot)/(4.0_dp*pi)
2819 END DO
2820 END DO
2821 ELSE IF (index(line, "LSD") /= 0) THEN
2822 ! LSD potential
2823 potential%lsdpot = .true.
2824 CALL parser_get_object(parser, potential%nexp_lsd)
2825 DO ipot = 1, potential%nexp_lsd
2826 CALL parser_get_object(parser, potential%alpha_lsd(ipot), newline=.true.)
2827 CALL parser_get_object(parser, potential%nct_lsd(ipot))
2828 DO ic = 1, potential%nct_lsd(ipot)
2829 CALL parser_get_object(parser, potential%cval_lsd(ic, ipot))
2830 END DO
2831 END DO
2832 ELSE
2833 cpabort("Parsing of extended potential type failed.")
2834 END IF
2835 ELSE
2836 cpabort("Invalid input token found.")
2837 END IF
2838 END DO
2839 ! Read the parameters for the non-local part of the GTH pseudopotential (ppnl)
2840 CALL parser_get_object(parser, nlmax)
2841 IF (nlmax > 0) THEN
2842 IF (parser_test_next_token(parser) == "STR") THEN
2843 CALL parser_get_object(parser, line)
2844 IF (index(line, "SOC") /= 0) potential%soc = .true.
2845 END IF
2846 ! Load the parameter for n non-local projectors
2847 DO l = 0, nlmax - 1
2848 CALL parser_get_object(parser, potential%rcnl(l), newline=.true.)
2849 CALL parser_get_object(parser, potential%nl(l))
2850 DO i = 1, potential%nl(l)
2851 IF (i == 1) THEN
2852 CALL parser_get_object(parser, potential%hnl(i, i, l))
2853 ELSE
2854 CALL parser_get_object(parser, potential%hnl(i, i, l), newline=.true.)
2855 END IF
2856 DO j = i + 1, potential%nl(l)
2857 CALL parser_get_object(parser, potential%hnl(i, j, l))
2858 potential%hnl(j, i, l) = potential%hnl(i, j, l)
2859 END DO
2860 END DO
2861 IF (potential%soc .AND. l /= 0) THEN
2862 DO i = 1, potential%nl(l)
2863 CALL parser_get_object(parser, potential%knl(i, i, l), newline=.true.)
2864 DO j = i + 1, potential%nl(l)
2865 CALL parser_get_object(parser, potential%knl(i, j, l))
2866 potential%knl(j, i, l) = potential%knl(i, j, l)
2867 END DO
2868 END DO
2869 END IF
2870 END DO
2871 END IF
2872 EXIT search_loop
2873 END IF
2874 ELSE
2875 ! Stop program, if the end of file is reached
2876 cpabort("End of file reached unexpectedly")
2877 END IF
2878
2879 END DO search_loop
2880
2881 CALL parser_release(parser)
2882 END block
2883 END IF
2884
2885 END SUBROUTINE read_gth_potential
2886! **************************************************************************************************
2887!> \brief ...
2888!> \param element_symbol ...
2889!> \param potential ...
2890!> \param pseudo_name ...
2891!> \param pseudo_file ...
2892!> \param potential_section ...
2893!> \param potential_found ...
2894! **************************************************************************************************
2895 SUBROUTINE read_ecp_potential_file(element_symbol, potential, pseudo_name, pseudo_file, &
2896 potential_section, potential_found)
2897
2898 CHARACTER(LEN=*), INTENT(IN) :: element_symbol
2899 TYPE(atom_ecppot_type), INTENT(INOUT) :: potential
2900 CHARACTER(LEN=*), INTENT(IN) :: pseudo_name, pseudo_file
2901 TYPE(section_vals_type), POINTER :: potential_section
2902 LOGICAL, INTENT(OUT), OPTIONAL :: potential_found
2903
2904 CHARACTER(LEN=240) :: line
2905 CHARACTER(len=5*default_string_length) :: line_att
2906 CHARACTER(LEN=LEN(element_symbol)+1) :: symbol
2907 CHARACTER(LEN=LEN(pseudo_name)) :: apname
2908 INTEGER :: i, ic, l, ncore, nel
2909 LOGICAL :: found, is_ok, read_from_input
2910 TYPE(cp_sll_val_type), POINTER :: list
2911 TYPE(val_type), POINTER :: val
2912
2913 apname = pseudo_name
2914 symbol = element_symbol
2915 IF (PRESENT(potential_found)) potential_found = .false.
2916 CALL get_ptable_info(symbol, number=ncore)
2917
2918 potential%symbol = symbol
2919 potential%pname = apname
2920 potential%econf = 0
2921 potential%zion = 0
2922 potential%lmax = -1
2923 potential%nloc = 0
2924 potential%nrloc = 0
2925 potential%aloc = 0.0_dp
2926 potential%bloc = 0.0_dp
2927 potential%npot = 0
2928 potential%nrpot = 0
2929 potential%apot = 0.0_dp
2930 potential%bpot = 0.0_dp
2931
2932 read_from_input = .false.
2933 CALL section_vals_get(potential_section, explicit=read_from_input)
2934 IF (read_from_input) THEN
2935 CALL section_vals_list_get(potential_section, "_DEFAULT_KEYWORD_", list=list)
2936 ! number of electrons (mandatory line)
2937 is_ok = cp_sll_val_next(list, val)
2938 cpassert(is_ok)
2939 CALL val_get(val, c_val=line_att)
2940 CALL remove_word(line_att)
2941 CALL remove_word(line_att)
2942 ! read number of electrons
2943 READ (line_att, *) nel
2944 potential%zion = real(ncore - nel, kind=dp)
2945 ! local potential (mandatory block)
2946 is_ok = cp_sll_val_next(list, val)
2947 cpassert(is_ok)
2948 CALL val_get(val, c_val=line_att)
2949 DO i = 1, 10
2950 IF (.NOT. cp_sll_val_next(list, val)) EXIT
2951 CALL val_get(val, c_val=line_att)
2952 IF (index(line_att, element_symbol) == 0) THEN
2953 potential%nloc = potential%nloc + 1
2954 ic = potential%nloc
2955 READ (line_att, *) potential%nrloc(ic), potential%bloc(ic), potential%aloc(ic)
2956 ELSE
2957 EXIT
2958 END IF
2959 END DO
2960 ! read potentials
2961 DO
2962 CALL val_get(val, c_val=line_att)
2963 IF (index(line_att, element_symbol) == 0) THEN
2964 potential%npot(l) = potential%npot(l) + 1
2965 ic = potential%npot(l)
2966 READ (line_att, *) potential%nrpot(ic, l), potential%bpot(ic, l), potential%apot(ic, l)
2967 ELSE
2968 potential%lmax = potential%lmax + 1
2969 l = potential%lmax
2970 END IF
2971 IF (.NOT. cp_sll_val_next(list, val)) EXIT
2972 END DO
2973
2974 ELSE
2975 block
2976 TYPE(cp_parser_type) :: parser
2977 CALL parser_create(parser, pseudo_file)
2978
2979 search_loop: DO
2980 CALL parser_search_string(parser, trim(apname), .true., found, line)
2981 IF (found) THEN
2982 match_loop: DO
2983 CALL parser_get_object(parser, line, newline=.true.)
2984 IF (trim(line) == element_symbol) THEN
2985 IF (PRESENT(potential_found)) potential_found = .true.
2986 CALL parser_get_object(parser, line, lower_to_upper=.true.)
2987 cpassert(trim(line) == "NELEC")
2988 ! read number of electrons
2989 CALL parser_get_object(parser, nel)
2990 potential%zion = real(ncore - nel, kind=dp)
2991 ! read local potential flag line "<XX> ul"
2992 CALL parser_get_object(parser, line, newline=.true.)
2993 ! read local potential
2994 DO i = 1, 15
2995 CALL parser_read_line(parser, 1)
2996 IF (parser_test_next_token(parser) == "STR") EXIT
2997 potential%nloc = potential%nloc + 1
2998 ic = potential%nloc
2999 CALL parser_get_object(parser, potential%nrloc(ic))
3000 CALL parser_get_object(parser, potential%bloc(ic))
3001 CALL parser_get_object(parser, potential%aloc(ic))
3002 END DO
3003 ! read potentials (start with l loop)
3004 DO l = 0, 15
3005 CALL parser_get_object(parser, symbol)
3006 IF (symbol == element_symbol) THEN
3007 ! new l block
3008 potential%lmax = potential%lmax + 1
3009 DO i = 1, 15
3010 CALL parser_read_line(parser, 1)
3011 IF (parser_test_next_token(parser) == "STR") EXIT
3012 potential%npot(l) = potential%npot(l) + 1
3013 ic = potential%npot(l)
3014 CALL parser_get_object(parser, potential%nrpot(ic, l))
3015 CALL parser_get_object(parser, potential%bpot(ic, l))
3016 CALL parser_get_object(parser, potential%apot(ic, l))
3017 END DO
3018 ELSE
3019 EXIT
3020 END IF
3021 END DO
3022 EXIT search_loop
3023 ELSE IF (line == "END") THEN
3024 IF (PRESENT(potential_found)) THEN
3025 CALL parser_release(parser)
3026 RETURN
3027 END IF
3028 cpabort("Element not found in ECP library")
3029 END IF
3030 END DO match_loop
3031 ELSE
3032 IF (PRESENT(potential_found)) THEN
3033 CALL parser_release(parser)
3034 RETURN
3035 END IF
3036 cpabort("ECP type not found in library")
3037 END IF
3038
3039 END DO search_loop
3040
3041 CALL parser_release(parser)
3042 END block
3043 END IF
3044
3045 IF (read_from_input .AND. PRESENT(potential_found)) potential_found = .true.
3046
3047 ! set up econf
3048 potential%econf(0:3) = ptable(ncore)%e_conv(0:3)
3049 SELECT CASE (nel)
3050 CASE DEFAULT
3051 cpabort("Unknown Core State")
3052 CASE (0)
3053 CASE (2)
3054 potential%econf(0:3) = potential%econf(0:3) - ptable(2)%e_conv(0:3)
3055 CASE (10)
3056 potential%econf(0:3) = potential%econf(0:3) - ptable(10)%e_conv(0:3)
3057 CASE (18)
3058 potential%econf(0:3) = potential%econf(0:3) - ptable(18)%e_conv(0:3)
3059 CASE (28)
3060 potential%econf(0:3) = potential%econf(0:3) - ptable(18)%e_conv(0:3)
3061 potential%econf(2) = potential%econf(2) - 10
3062 CASE (36)
3063 potential%econf(0:3) = potential%econf(0:3) - ptable(36)%e_conv(0:3)
3064 CASE (46)
3065 potential%econf(0:3) = potential%econf(0:3) - ptable(36)%e_conv(0:3)
3066 potential%econf(2) = potential%econf(2) - 10
3067 CASE (54)
3068 potential%econf(0:3) = potential%econf(0:3) - ptable(54)%e_conv(0:3)
3069 CASE (60)
3070 potential%econf(0:3) = potential%econf(0:3) - ptable(36)%e_conv(0:3)
3071 potential%econf(2) = potential%econf(2) - 10
3072 potential%econf(3) = potential%econf(3) - 14
3073 CASE (68)
3074 potential%econf(0:3) = potential%econf(0:3) - ptable(54)%e_conv(0:3)
3075 potential%econf(3) = potential%econf(3) - 14
3076 CASE (78)
3077 potential%econf(0:3) = potential%econf(0:3) - ptable(54)%e_conv(0:3)
3078 potential%econf(2) = potential%econf(2) - 10
3079 potential%econf(3) = potential%econf(3) - 14
3080 END SELECT
3081 !
3082 cpassert(all(potential%econf >= 0))
3083
3084 END SUBROUTINE read_ecp_potential_file
3085
3086! **************************************************************************************************
3087!> \brief Read an ECP potential from an ordered list of library files.
3088!> \param element_symbol ...
3089!> \param potential ...
3090!> \param pseudo_name ...
3091!> \param pseudo_files ...
3092!> \param potential_section ...
3093! **************************************************************************************************
3094 SUBROUTINE read_ecp_potential_files(element_symbol, potential, pseudo_name, pseudo_files, &
3095 potential_section)
3096
3097 CHARACTER(LEN=*), INTENT(IN) :: element_symbol
3098 TYPE(atom_ecppot_type), INTENT(INOUT) :: potential
3099 CHARACTER(LEN=*), INTENT(IN) :: pseudo_name
3100 CHARACTER(LEN=*), DIMENSION(:), INTENT(IN) :: pseudo_files
3101 TYPE(section_vals_type), POINTER :: potential_section
3102
3103 CHARACTER(LEN=:), ALLOCATABLE :: file_list
3104 INTEGER :: i
3105 LOGICAL :: potential_found
3106
3107 DO i = 1, SIZE(pseudo_files)
3108 CALL read_ecp_potential_file(element_symbol, potential, pseudo_name, pseudo_files(i), &
3109 potential_section, potential_found)
3110 IF (potential_found) RETURN
3111 END DO
3112
3113 file_list = ""
3114 DO i = 1, SIZE(pseudo_files)
3115 file_list = trim(file_list)//"<"//trim(pseudo_files(i))//"> "
3116 END DO
3117 CALL cp_abort(__location__, &
3118 "The requested ECP potential <"//trim(pseudo_name)// &
3119 "> for element <"//trim(element_symbol)// &
3120 "> was not found in the potential files "//trim(file_list))
3121
3122 END SUBROUTINE read_ecp_potential_files
3123! **************************************************************************************************
3124!> \brief ...
3125!> \param grid1 ...
3126!> \param grid2 ...
3127!> \return ...
3128! **************************************************************************************************
3129 FUNCTION atom_compare_grids(grid1, grid2) RESULT(is_equal)
3130 TYPE(grid_atom_type) :: grid1, grid2
3131 LOGICAL :: is_equal
3132
3133 INTEGER :: i
3134 REAL(kind=dp) :: dr, dw
3135
3136 is_equal = .true.
3137 IF (grid1%nr == grid2%nr) THEN
3138 DO i = 1, grid2%nr
3139 dr = abs(grid1%rad(i) - grid2%rad(i))
3140 dw = abs(grid1%wr(i) - grid2%wr(i))
3141 IF (dr + dw > 1.0e-12_dp) THEN
3142 is_equal = .false.
3143 EXIT
3144 END IF
3145 END DO
3146 ELSE
3147 is_equal = .false.
3148 END IF
3149
3150 END FUNCTION atom_compare_grids
3151! **************************************************************************************************
3152
3153END MODULE atom_types
Define the atom type and its sub types.
Definition atom_types.F:15
integer, parameter, public num_basis
Definition atom_types.F:70
subroutine, public read_atom_opt_section(optimization, opt_section)
...
subroutine, public create_atom_type(atom)
...
Definition atom_types.F:968
integer, parameter, public cgto_basis
Definition atom_types.F:70
integer, parameter, public gto_basis
Definition atom_types.F:70
integer, parameter, public sto_basis
Definition atom_types.F:70
subroutine, public release_atom_type(atom)
...
Definition atom_types.F:992
subroutine, public release_opmat(opmat)
...
subroutine, public release_atom_potential(potential)
...
subroutine, public init_atom_basis(basis, basis_section, zval, btyp)
Initialize the basis for the atomic code.
Definition atom_types.F:384
logical function, public atom_compare_grids(grid1, grid2)
...
subroutine, public release_opgrid(opgrid)
...
subroutine, public release_atom_orbs(orbs)
...
integer, parameter, public lmat
Definition atom_types.F:68
subroutine, public set_atom(atom, basis, state, integrals, orbitals, potential, zcore, pp_calc, do_zmp, doread, read_vxc, method_type, relativistic, coulomb_integral_type, exchange_integral_type, fmat)
...
subroutine, public init_atom_potential(potential, potential_section, zval)
...
subroutine, public release_atom_basis(basis)
...
Definition atom_types.F:934
subroutine, public create_atom_orbs(orbs, mbas, mo)
...
subroutine, public init_atom_basis_default_pp(basis)
...
Definition atom_types.F:723
subroutine, public create_opmat(opmat, n, lmax)
...
subroutine, public atom_basis_gridrep(basis, gbasis, r, rab)
...
Definition atom_types.F:800
subroutine, public setup_hf_section(hf_frac, do_hfx, atom, xc_section, extype)
...
subroutine, public clementi_geobas(zval, cval, aval, ngto, ival)
...
subroutine, public create_opgrid(opgrid, grid)
...
Routines that process Quantum Espresso UPF files.
Definition atom_upf.F:14
pure subroutine, public atom_release_upf(upfpot)
...
Definition atom_upf.F:1191
subroutine, public atom_read_upf(pot, upf_filename, read_header, read_augmentation)
...
Definition atom_upf.F:52
Definition atom.F:9
Calculates Bessel functions.
Definition bessel_lib.F:16
elemental impure real(kind=dp) function, public bessel0(x, l)
...
Definition bessel_lib.F:161
collects all references to literature in CP2K as new algorithms / method are included from literature...
integer, save, public limpanuparb2011
logical function, public cp_sll_val_next(iterator, el_att)
returns true if the actual element is valid (i.e. iterator ont at end) moves the iterator to the next...
Utility routines to read data from files. Kept as close as possible to the old parser because.
subroutine, public parser_read_line(parser, nline, at_end)
Read the next line from a logical unit "unit" (I/O node only). Skip (nline-1) lines and skip also all...
subroutine, public parser_get_next_line(parser, nline, at_end)
Read the next input line and broadcast the input information. Skip (nline-1) lines and skip also all ...
character(len=3) function, public parser_test_next_token(parser, string_length)
Test next input object.
subroutine, public parser_search_string(parser, string, ignore_case, found, line, begin_line, search_from_begin_of_file)
Search a string pattern in a file defined by its logical unit number "unit". A case sensitive search ...
Utility routines to read data from files. Kept as close as possible to the old parser because.
subroutine, public parser_release(parser)
releases the parser
subroutine, public parser_create(parser, file_name, unit_nr, para_env, end_section_label, separator_chars, comment_char, continuation_char, quote_char, section_char, parse_white_lines, initial_variables, apply_preprocessing)
Start a parser run. Initial variables allow to @SET stuff before opening the file.
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public do_rks_atom
integer, parameter, public do_analytic
integer, parameter, public sgp_pseudo
integer, parameter, public gth_pseudo
integer, parameter, public ecp_pseudo
integer, parameter, public do_nonrel_atom
integer, parameter, public no_conf
integer, parameter, public do_potential_mix_cl
integer, parameter, public upf_pseudo
integer, parameter, public contracted_gto
integer, parameter, public poly_conf
integer, parameter, public no_pseudo
integer, parameter, public barrier_conf
integer, parameter, public do_numeric
integer, parameter, public do_gapw_log
integer, parameter, public do_potential_coulomb
integer, parameter, public gaussian
integer, parameter, public do_potential_short
integer, parameter, public do_semi_analytic
integer, parameter, public geometrical_gto
integer, parameter, public do_potential_long
integer, parameter, public numerical
integer, parameter, public slater
objects that represent the structure of input sections and the data contained in an input section
subroutine, public section_vals_list_get(section_vals, keyword_name, i_rep_section, list)
returns the requested list
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
a wrapper for basic fortran types.
subroutine, public val_get(val, has_l, has_i, has_r, has_lc, has_c, l_val, l_vals, i_val, i_vals, r_val, r_vals, c_val, c_vals, len_c, type_of_var, enum)
returns the stored values
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
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition list.F:24
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
real(kind=dp), dimension(-1:2 *maxfac+1), parameter, public dfac
real(kind=dp), parameter, public rootpi
real(kind=dp), dimension(0:maxfac), parameter, public fac
Periodic Table related data definitions.
type(atom), dimension(0:nelem), public ptable
subroutine, public get_ptable_info(symbol, number, amass, ielement, covalent_radius, metallic_radius, vdw_radius, found)
Pass information about the kind given the element symbol.
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 create_grid_atom(grid_atom, nr, na, llmax, ll, quadrature)
...
Utilities for string manipulations.
character(len=1), parameter, public newline
subroutine, public remove_word(string)
remove a word from a string (words are separated by white spaces)
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.
Provides all information about a basis set.
Definition atom_types.F:79
Provides all information about a pseudopotential.
Definition atom_types.F:99
Provides info about hartree-fock exchange (For now, we only support potentials that can be represente...
Definition atom_types.F:188
Information on optimization procedure.
Definition atom_types.F:284
Holds atomic orbitals and energies.
Definition atom_types.F:238
Provides all information on states and occupation.
Definition atom_types.F:199
Provides all information about an atomic kind.
Definition atom_types.F:294
Holds atomic integrals.
Definition atom_types.F:213
Operator grids.
Definition atom_types.F:259
Operator matrices.
Definition atom_types.F:252
represent a single linked list that stores pointers to the elements
a type to have a wrapper that stores any basic fortran type