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