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