(git:f56c6e3)
Loading...
Searching...
No Matches
atom_kind_orbitals.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief calculate the orbitals for a given atomic kind type
10! **************************************************************************************************
11MODULE atom_kind_orbitals
12 USE ai_onecenter, ONLY: sg_erfc
13 USE atom_electronic_structure, ONLY: calculate_atom
14 USE atom_fit, ONLY: atom_fit_density
15 USE atom_operators, ONLY: atom_int_release,&
16 atom_int_setup,&
17 atom_ppint_release,&
18 atom_ppint_setup,&
19 atom_relint_release,&
20 atom_relint_setup
21 USE atom_set_basis, ONLY: set_kind_basis_atomic
22 USE atom_types, ONLY: &
23 cgto_basis, clementi_geobas, gto_basis, atom_basis_type, atom_ecppot_type, &
24 atom_gthpot_type, atom_integrals, atom_orbitals, atom_potential_type, atom_sgppot_type, &
25 atom_type, create_atom_orbs, create_atom_type, lmat, release_atom_basis, &
26 release_atom_potential, release_atom_type, set_atom
27 USE atom_utils, ONLY: atom_density,&
28 get_maxl_occ,&
29 get_maxn_occ
30 USE atomic_kind_types, ONLY: atomic_kind_type,&
31 get_atomic_kind
32 USE basis_set_types, ONLY: get_gto_basis_set,&
33 gto_basis_set_type
34 USE external_potential_types, ONLY: all_potential_type,&
35 get_potential,&
36 gth_potential_type,&
37 sgp_potential_type
38 USE input_constants, ONLY: &
39 barrier_conf, do_analytic, do_dkh0_atom, do_dkh1_atom, do_dkh2_atom, do_dkh3_atom, &
40 do_gapw_log, do_nonrel_atom, do_numeric, do_rks_atom, do_sczoramp_atom, do_uks_atom, &
41 do_zoramp_atom, ecp_pseudo, gth_pseudo, no_pseudo, poly_conf, rel_dkh, rel_none, &
42 rel_sczora_mp, rel_zora, rel_zora_full, rel_zora_mp, sgp_pseudo
43 USE input_section_types, ONLY: section_vals_type
44 USE kinds, ONLY: dp
45 USE mathconstants, ONLY: dfac,&
46 pi
47 USE periodic_table, ONLY: ptable
48 USE physcon, ONLY: bohr
49 USE qs_grid_atom, ONLY: allocate_grid_atom,&
50 create_grid_atom,&
51 grid_atom_type
52 USE qs_kind_types, ONLY: get_qs_kind,&
53 init_atom_electronic_state,&
54 qs_kind_type,&
55 set_pseudo_state
56 USE rel_control_types, ONLY: rel_control_type
57#include "./base/base_uses.f90"
58
59 IMPLICIT NONE
60
61 PRIVATE
62
63 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_kind_orbitals'
64
65 PUBLIC :: calculate_atomic_orbitals, calculate_atomic_density, &
66 calculate_atomic_relkin, gth_potential_conversion
67
68! **************************************************************************************************
69
70CONTAINS
71
72! **************************************************************************************************
73!> \brief ...
74!> \param atomic_kind ...
75!> \param qs_kind ...
76!> \param agrid ...
77!> \param iunit ...
78!> \param pmat ...
79!> \param fmat ...
80!> \param density ...
81!> \param wavefunction ...
82!> \param wfninfo ...
83!> \param confine ...
84!> \param xc_section ...
85!> \param nocc ...
86! **************************************************************************************************
87 SUBROUTINE calculate_atomic_orbitals(atomic_kind, qs_kind, agrid, iunit, pmat, fmat, &
88 density, wavefunction, wfninfo, confine, xc_section, nocc)
89 TYPE(atomic_kind_type), INTENT(IN) :: atomic_kind
90 TYPE(qs_kind_type), INTENT(IN) :: qs_kind
91 TYPE(grid_atom_type), OPTIONAL :: agrid
92 INTEGER, INTENT(IN), OPTIONAL :: iunit
93 REAL(kind=dp), DIMENSION(:, :, :), OPTIONAL, &
94 POINTER :: pmat, fmat
95 REAL(kind=dp), DIMENSION(:), OPTIONAL, POINTER :: density
96 REAL(kind=dp), DIMENSION(:, :), OPTIONAL, POINTER :: wavefunction, wfninfo
97 LOGICAL, INTENT(IN), OPTIONAL :: confine
98 TYPE(section_vals_type), OPTIONAL, POINTER :: xc_section
99 INTEGER, DIMENSION(:), OPTIONAL :: nocc
100
101 INTEGER :: i, ii, j, k, k1, k2, l, ll, m, mb, mo, &
102 nr, nset, nsgf, z
103 INTEGER, DIMENSION(0:lmat) :: nbb
104 INTEGER, DIMENSION(0:lmat, 10) :: ncalc, ncore, nelem
105 INTEGER, DIMENSION(0:lmat, 100) :: set_index, shell_index
106 INTEGER, DIMENSION(:), POINTER :: nshell
107 INTEGER, DIMENSION(:, :), POINTER :: first_sgf, ls
108 LOGICAL :: ecp_semi_local, ghost, has_pp, uks
109 REAL(kind=dp) :: ok, scal, zeff
110 REAL(kind=dp), DIMENSION(0:lmat, 10, 2) :: edelta
111 TYPE(all_potential_type), POINTER :: all_potential
112 TYPE(atom_basis_type), POINTER :: basis
113 TYPE(atom_integrals), POINTER :: integrals
114 TYPE(atom_orbitals), POINTER :: orbitals
115 TYPE(atom_potential_type), POINTER :: potential
116 TYPE(atom_type), POINTER :: atom
117 TYPE(gth_potential_type), POINTER :: gth_potential
118 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
119 TYPE(sgp_potential_type), POINTER :: sgp_potential
120
121 NULLIFY (atom)
122 CALL create_atom_type(atom)
123
124 IF (PRESENT(xc_section)) THEN
125 atom%xc_section => xc_section
126 ELSE
127 NULLIFY (atom%xc_section)
128 END IF
129
130 CALL get_atomic_kind(atomic_kind, z=z)
131 NULLIFY (all_potential, gth_potential, sgp_potential, orb_basis_set)
132 CALL get_qs_kind(qs_kind, zeff=zeff, &
133 basis_set=orb_basis_set, &
134 ghost=ghost, &
135 all_potential=all_potential, &
136 gth_potential=gth_potential, &
137 sgp_potential=sgp_potential)
138
139 has_pp = ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)
140
141 atom%z = z
142 CALL set_atom(atom, &
143 pp_calc=has_pp, &
144 do_zmp=.false., &
145 doread=.false., &
146 read_vxc=.false., &
147 relativistic=do_nonrel_atom, &
148 coulomb_integral_type=do_numeric, &
149 exchange_integral_type=do_numeric)
150
151 ALLOCATE (potential, integrals)
152
153 IF (PRESENT(confine)) THEN
154 potential%confinement = confine
155 ELSE
156 IF (ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)) THEN
157 potential%confinement = .true.
158 ELSE
159 potential%confinement = .false.
160 END IF
161 END IF
162 potential%conf_type = poly_conf
163 potential%acon = 0.1_dp
164 potential%rcon = 2.0_dp*ptable(z)%vdw_radius*bohr
165 potential%scon = 2.0_dp
166
167 IF (ASSOCIATED(gth_potential)) THEN
168 potential%ppot_type = gth_pseudo
169 CALL get_potential(gth_potential, zeff=zeff)
170 CALL gth_potential_conversion(gth_potential, potential%gth_pot)
171 CALL set_atom(atom, zcore=nint(zeff), potential=potential)
172 ELSE IF (ASSOCIATED(sgp_potential)) THEN
173 CALL get_potential(sgp_potential, ecp_semi_local=ecp_semi_local)
174 IF (ecp_semi_local) THEN
175 potential%ppot_type = ecp_pseudo
176 CALL ecp_potential_conversion(sgp_potential, potential%ecp_pot)
177 potential%ecp_pot%symbol = ptable(z)%symbol
178 ELSE
179 potential%ppot_type = sgp_pseudo
180 CALL sgp_potential_conversion(sgp_potential, potential%sgp_pot)
181 potential%sgp_pot%symbol = ptable(z)%symbol
182 END IF
183 CALL get_potential(sgp_potential, zeff=zeff)
184 CALL set_atom(atom, zcore=nint(zeff), potential=potential)
185 ELSE
186 potential%ppot_type = no_pseudo
187 CALL set_atom(atom, zcore=z, potential=potential)
188 END IF
189
190 NULLIFY (basis)
191 ALLOCATE (basis)
192 CALL set_kind_basis_atomic(basis, orb_basis_set, has_pp, agrid)
193
194 CALL set_atom(atom, basis=basis)
195
196 ! optimization defaults
197 atom%optimization%damping = 0.2_dp
198 atom%optimization%eps_scf = 1.e-6_dp
199 atom%optimization%eps_diis = 100._dp
200 atom%optimization%max_iter = 50
201 atom%optimization%n_diis = 5
202
203 ! set up the electronic state
204 CALL init_atom_electronic_state(atomic_kind=atomic_kind, &
205 qs_kind=qs_kind, &
206 ncalc=ncalc, &
207 ncore=ncore, &
208 nelem=nelem, &
209 edelta=edelta)
210
211 ! restricted or unrestricted?
212 IF (sum(abs(edelta)) > 0.0_dp) THEN
213 uks = .true.
214 CALL set_atom(atom, method_type=do_uks_atom)
215 ELSE
216 uks = .false.
217 CALL set_atom(atom, method_type=do_rks_atom)
218 END IF
219
220 ALLOCATE (atom%state)
221
222 atom%state%core = 0._dp
223 atom%state%core(0:lmat, 1:7) = real(ncore(0:lmat, 1:7), dp)
224 atom%state%occ = 0._dp
225 IF (uks) THEN
226 atom%state%occ(0:lmat, 1:7) = real(ncalc(0:lmat, 1:7), dp) + &
227 edelta(0:lmat, 1:7, 1) + edelta(0:lmat, 1:7, 2)
228 ELSE
229 atom%state%occ(0:lmat, 1:7) = real(ncalc(0:lmat, 1:7), dp)
230 END IF
231 atom%state%occupation = 0._dp
232 DO l = 0, lmat
233 k = 0
234 DO i = 1, 7
235 IF (ncalc(l, i) > 0) THEN
236 k = k + 1
237 IF (uks) THEN
238 atom%state%occupation(l, k) = real(ncalc(l, i), dp) + &
239 edelta(l, i, 1) + edelta(l, i, 2)
240 atom%state%occa(l, k) = 0.5_dp*real(ncalc(l, i), dp) + edelta(l, i, 1)
241 atom%state%occb(l, k) = 0.5_dp*real(ncalc(l, i), dp) + edelta(l, i, 2)
242 ELSE
243 atom%state%occupation(l, k) = real(ncalc(l, i), dp)
244 END IF
245 END IF
246 END DO
247 ok = real(2*l + 1, kind=dp)
248 IF (uks) THEN
249 DO i = 1, 7
250 atom%state%occ(l, i) = min(atom%state%occ(l, i), 2.0_dp*ok)
251 atom%state%occa(l, i) = min(atom%state%occa(l, i), ok)
252 atom%state%occb(l, i) = min(atom%state%occb(l, i), ok)
253 atom%state%occupation(l, i) = atom%state%occa(l, i) + atom%state%occb(l, i)
254 END DO
255 ELSE
256 DO i = 1, 7
257 atom%state%occ(l, i) = min(atom%state%occ(l, i), 2.0_dp*ok)
258 atom%state%occupation(l, i) = min(atom%state%occupation(l, i), 2.0_dp*ok)
259 END DO
260 END IF
261 END DO
262 IF (uks) THEN
263 atom%state%multiplicity = nint(abs(sum(atom%state%occa - atom%state%occb)) + 1)
264 ELSE
265 atom%state%multiplicity = -1
266 END IF
267
268 atom%state%maxl_occ = get_maxl_occ(atom%state%occupation)
269 atom%state%maxn_occ = get_maxn_occ(atom%state%occupation)
270 atom%state%maxl_calc = atom%state%maxl_occ
271 atom%state%maxn_calc = atom%state%maxn_occ
272
273 ! total number of occupied orbitals
274 IF (PRESENT(nocc) .AND. ghost) THEN
275 nocc = 0
276 ELSEIF (PRESENT(nocc)) THEN
277 nocc = 0
278 DO l = 0, lmat
279 DO k = 1, 7
280 IF (uks) THEN
281 IF (atom%state%occa(l, k) > 0.0_dp) THEN
282 nocc(1) = nocc(1) + 2*l + 1
283 END IF
284 IF (atom%state%occb(l, k) > 0.0_dp) THEN
285 nocc(2) = nocc(2) + 2*l + 1
286 END IF
287 ELSE
288 IF (atom%state%occupation(l, k) > 0.0_dp) THEN
289 nocc(1) = nocc(1) + 2*l + 1
290 nocc(2) = nocc(2) + 2*l + 1
291 END IF
292 END IF
293 END DO
294 END DO
295 END IF
296
297 ! calculate integrals
298 ! general integrals
299 CALL atom_int_setup(integrals, basis, potential=atom%potential, &
300 eri_coulomb=(atom%coulomb_integral_type == do_analytic), &
301 eri_exchange=(atom%exchange_integral_type == do_analytic))
302 ! potential
303 CALL atom_ppint_setup(integrals, basis, potential=atom%potential)
304 ! relativistic correction terms
305 NULLIFY (integrals%tzora, integrals%hdkh)
306 CALL atom_relint_setup(integrals, basis, atom%relativistic, zcore=real(atom%zcore, dp))
307 CALL set_atom(atom, integrals=integrals)
308
309 NULLIFY (orbitals)
310 mo = maxval(atom%state%maxn_calc)
311 mb = maxval(atom%basis%nbas)
312 CALL create_atom_orbs(orbitals, mb, mo)
313 CALL set_atom(atom, orbitals=orbitals)
314
315 IF (.NOT. ghost) THEN
316 IF (PRESENT(iunit)) THEN
317 CALL calculate_atom(atom, iunit)
318 ELSE
319 CALL calculate_atom(atom, -1)
320 END IF
321 END IF
322
323 IF (PRESENT(pmat)) THEN
324 ! recover density matrix in CP2K/GPW order and normalization
325 CALL get_gto_basis_set(orb_basis_set, &
326 nset=nset, nshell=nshell, l=ls, nsgf=nsgf, first_sgf=first_sgf)
327 set_index = 0
328 shell_index = 0
329 nbb = 0
330 DO i = 1, nset
331 DO j = 1, nshell(i)
332 l = ls(j, i)
333 IF (l <= lmat) THEN
334 nbb(l) = nbb(l) + 1
335 k = nbb(l)
336 cpassert(k <= 100)
337 set_index(l, k) = i
338 shell_index(l, k) = j
339 END IF
340 END DO
341 END DO
342
343 IF (ASSOCIATED(pmat)) THEN
344 DEALLOCATE (pmat)
345 END IF
346 ALLOCATE (pmat(nsgf, nsgf, 2))
347 pmat = 0._dp
348 IF (.NOT. ghost) THEN
349 DO l = 0, lmat
350 ll = 2*l
351 DO k1 = 1, atom%basis%nbas(l)
352 DO k2 = 1, atom%basis%nbas(l)
353 scal = sqrt(atom%integrals%ovlp(k1, k1, l)*atom%integrals%ovlp(k2, k2, l))/real(2*l + 1, kind=dp)
354 i = first_sgf(shell_index(l, k1), set_index(l, k1))
355 j = first_sgf(shell_index(l, k2), set_index(l, k2))
356 IF (uks) THEN
357 DO m = 0, ll
358 pmat(i + m, j + m, 1) = atom%orbitals%pmata(k1, k2, l)*scal
359 pmat(i + m, j + m, 2) = atom%orbitals%pmatb(k1, k2, l)*scal
360 END DO
361 ELSE
362 DO m = 0, ll
363 pmat(i + m, j + m, 1) = atom%orbitals%pmat(k1, k2, l)*scal
364 END DO
365 END IF
366 END DO
367 END DO
368 END DO
369 IF (uks) THEN
370 pmat(:, :, 1) = pmat(:, :, 1) + pmat(:, :, 2)
371 pmat(:, :, 2) = pmat(:, :, 1) - 2.0_dp*pmat(:, :, 2)
372 END IF
373 END IF
374 END IF
375
376 IF (PRESENT(fmat)) THEN
377 ! recover fock matrix in CP2K/GPW order.
378 ! Caution: Normalization is not take care of, so it's probably weird.
379 CALL get_gto_basis_set(orb_basis_set, &
380 nset=nset, nshell=nshell, l=ls, nsgf=nsgf, first_sgf=first_sgf)
381 set_index = 0
382 shell_index = 0
383 nbb = 0
384 DO i = 1, nset
385 DO j = 1, nshell(i)
386 l = ls(j, i)
387 IF (l <= lmat) THEN
388 nbb(l) = nbb(l) + 1
389 k = nbb(l)
390 cpassert(k <= 100)
391 set_index(l, k) = i
392 shell_index(l, k) = j
393 END IF
394 END DO
395 END DO
396 IF (uks) cpabort("calculate_atomic_orbitals: only RKS is implemented")
397 IF (ASSOCIATED(fmat)) cpabort("fmat already associated")
398 IF (.NOT. ASSOCIATED(atom%fmat)) cpabort("atom%fmat not associated")
399 ALLOCATE (fmat(nsgf, nsgf, 1))
400 fmat = 0.0_dp
401 IF (.NOT. ghost) THEN
402 DO l = 0, lmat
403 ll = 2*l
404 DO k1 = 1, atom%basis%nbas(l)
405 DO k2 = 1, atom%basis%nbas(l)
406 scal = sqrt(atom%integrals%ovlp(k1, k1, l)*atom%integrals%ovlp(k2, k2, l))
407 i = first_sgf(shell_index(l, k1), set_index(l, k1))
408 j = first_sgf(shell_index(l, k2), set_index(l, k2))
409 DO m = 0, ll
410 fmat(i + m, j + m, 1) = atom%fmat%op(k1, k2, l)/scal
411 END DO
412 END DO
413 END DO
414 END DO
415 END IF
416 END IF
417
418 nr = basis%grid%nr
419
420 IF (PRESENT(density)) THEN
421 IF (ASSOCIATED(density)) DEALLOCATE (density)
422 ALLOCATE (density(nr))
423 IF (ghost) THEN
424 density = 0.0_dp
425 ELSE
426 CALL atom_density(density, atom%orbitals%pmat, atom%basis, atom%state%maxl_occ)
427 END IF
428 END IF
429
430 IF (PRESENT(wavefunction)) THEN
431 cpassert(PRESENT(wfninfo))
432 IF (ASSOCIATED(wavefunction)) DEALLOCATE (wavefunction)
433 IF (ASSOCIATED(wfninfo)) DEALLOCATE (wfninfo)
434 mo = sum(atom%state%maxn_occ)
435 ALLOCATE (wavefunction(nr, mo), wfninfo(2, mo))
436 wavefunction = 0.0_dp
437 IF (.NOT. ghost) THEN
438 ii = 0
439 DO l = 0, lmat
440 DO i = 1, atom%state%maxn_occ(l)
441 IF (atom%state%occupation(l, i) > 0.0_dp) THEN
442 ii = ii + 1
443 wfninfo(1, ii) = atom%state%occupation(l, i)
444 wfninfo(2, ii) = real(l, dp)
445 DO j = 1, atom%basis%nbas(l)
446 wavefunction(:, ii) = wavefunction(:, ii) + &
447 atom%orbitals%wfn(j, i, l)*basis%bf(:, j, l)
448 END DO
449 END IF
450 END DO
451 END DO
452 cpassert(mo == ii)
453 END IF
454 END IF
455
456 ! clean up
457 CALL atom_int_release(integrals)
458 CALL atom_ppint_release(integrals)
459 CALL atom_relint_release(integrals)
460 CALL release_atom_basis(basis)
461 CALL release_atom_potential(potential)
462 CALL release_atom_type(atom)
463
464 DEALLOCATE (potential, basis, integrals)
465
466 END SUBROUTINE calculate_atomic_orbitals
467
468! **************************************************************************************************
469!> \brief ...
470!> \param density ...
471!> \param atomic_kind ...
472!> \param qs_kind ...
473!> \param ngto ...
474!> \param iunit ...
475!> \param optbasis ... Default=T, if basis should be optimized, if not basis is given in input (density)
476!> \param allelectron ...
477!> \param confine ...
478! **************************************************************************************************
479 SUBROUTINE calculate_atomic_density(density, atomic_kind, qs_kind, ngto, iunit, &
480 optbasis, allelectron, confine)
481 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: density
482 TYPE(atomic_kind_type), POINTER :: atomic_kind
483 TYPE(qs_kind_type), POINTER :: qs_kind
484 INTEGER, INTENT(IN) :: ngto
485 INTEGER, INTENT(IN), OPTIONAL :: iunit
486 LOGICAL, INTENT(IN), OPTIONAL :: optbasis, allelectron, confine
487
488 INTEGER, PARAMETER :: num_gto = 40
489
490 INTEGER :: i, ii, iw, k, l, ll, m, mb, mo, ngp, nn, &
491 nr, quadtype, relativistic, z
492 INTEGER, DIMENSION(0:lmat) :: starti
493 INTEGER, DIMENSION(0:lmat, 10) :: ncalc, ncore, nelem
494 INTEGER, DIMENSION(:), POINTER :: econf
495 LOGICAL :: do_basopt, ecp_semi_local, monovalent
496 REAL(kind=dp) :: al, aval, cc, cval, ear, rk, xx, zeff
497 REAL(kind=dp), DIMENSION(num_gto+2) :: results
498 TYPE(all_potential_type), POINTER :: all_potential
499 TYPE(atom_basis_type), POINTER :: basis
500 TYPE(atom_integrals), POINTER :: integrals
501 TYPE(atom_orbitals), POINTER :: orbitals
502 TYPE(atom_potential_type), POINTER :: potential
503 TYPE(atom_type), POINTER :: atom
504 TYPE(grid_atom_type), POINTER :: grid
505 TYPE(gth_potential_type), POINTER :: gth_potential
506 TYPE(sgp_potential_type), POINTER :: sgp_potential
507
508 NULLIFY (atom)
509 CALL create_atom_type(atom)
510
511 CALL get_atomic_kind(atomic_kind, z=z)
512 NULLIFY (all_potential, gth_potential)
513 CALL get_qs_kind(qs_kind, zeff=zeff, &
514 all_potential=all_potential, &
515 gth_potential=gth_potential, &
516 sgp_potential=sgp_potential, &
517 monovalent=monovalent)
518
519 IF (PRESENT(iunit)) THEN
520 iw = iunit
521 ELSE
522 iw = -1
523 END IF
524
525 IF (PRESENT(allelectron)) THEN
526 IF (allelectron) THEN
527 NULLIFY (gth_potential)
528 zeff = z
529 END IF
530 END IF
531
532 do_basopt = .true.
533 IF (PRESENT(optbasis)) THEN
534 do_basopt = optbasis
535 END IF
536
537 cpassert(ngto <= num_gto)
538
539 IF (ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)) THEN
540 ! PP calculation are non-relativistic
541 relativistic = do_nonrel_atom
542 ELSE
543 ! AE calculations use DKH2
544 relativistic = do_dkh2_atom
545 END IF
546
547 atom%z = z
548 CALL set_atom(atom, &
549 pp_calc=(ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)), &
550 method_type=do_rks_atom, &
551 relativistic=relativistic, &
552 coulomb_integral_type=do_numeric, &
553 exchange_integral_type=do_numeric)
554
555 ALLOCATE (potential, basis, integrals)
556
557 IF (PRESENT(confine)) THEN
558 potential%confinement = confine
559 ELSE
560 IF (ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)) THEN
561 potential%confinement = .true.
562 ELSE
563 potential%confinement = .false.
564 END IF
565 END IF
566 potential%conf_type = barrier_conf
567 potential%acon = 200._dp
568 potential%rcon = 4.0_dp
569 potential%scon = 8.0_dp
570
571 IF (ASSOCIATED(gth_potential)) THEN
572 potential%ppot_type = gth_pseudo
573 CALL get_potential(gth_potential, zeff=zeff)
574 CALL gth_potential_conversion(gth_potential, potential%gth_pot)
575 CALL set_atom(atom, zcore=nint(zeff), potential=potential)
576 ELSE IF (ASSOCIATED(sgp_potential)) THEN
577 CALL get_potential(sgp_potential, ecp_semi_local=ecp_semi_local)
578 IF (ecp_semi_local) THEN
579 potential%ppot_type = ecp_pseudo
580 CALL ecp_potential_conversion(sgp_potential, potential%ecp_pot)
581 potential%ecp_pot%symbol = ptable(z)%symbol
582 ELSE
583 potential%ppot_type = sgp_pseudo
584 CALL sgp_potential_conversion(sgp_potential, potential%sgp_pot)
585 potential%sgp_pot%symbol = ptable(z)%symbol
586 END IF
587 CALL get_potential(sgp_potential, zeff=zeff)
588 CALL set_atom(atom, zcore=nint(zeff), potential=potential)
589 ELSE
590 potential%ppot_type = no_pseudo
591 CALL set_atom(atom, zcore=z, potential=potential)
592 END IF
593
594 ! atomic grid
595 NULLIFY (grid)
596 ngp = 400
597 quadtype = do_gapw_log
598 CALL allocate_grid_atom(grid)
599 CALL create_grid_atom(grid, ngp, 1, 1, 0, quadtype)
600 grid%nr = ngp
601 basis%grid => grid
602
603 NULLIFY (basis%am, basis%cm, basis%as, basis%ns, basis%bf, basis%dbf, basis%ddbf)
604
605 ! fill in the basis data structures
606 basis%eps_eig = 1.e-12_dp
607 basis%basis_type = gto_basis
608 CALL clementi_geobas(z, cval, aval, basis%nbas, starti)
609 basis%nprim = basis%nbas
610 m = maxval(basis%nbas)
611 ALLOCATE (basis%am(m, 0:lmat))
612 basis%am = 0._dp
613 DO l = 0, lmat
614 DO i = 1, basis%nbas(l)
615 ll = i - 1 + starti(l)
616 basis%am(i, l) = aval*cval**(ll)
617 END DO
618 END DO
619
620 basis%geometrical = .true.
621 basis%aval = aval
622 basis%cval = cval
623 basis%start = starti
624
625 ! initialize basis function on a radial grid
626 nr = basis%grid%nr
627 m = maxval(basis%nbas)
628 ALLOCATE (basis%bf(nr, m, 0:lmat))
629 ALLOCATE (basis%dbf(nr, m, 0:lmat))
630 ALLOCATE (basis%ddbf(nr, m, 0:lmat))
631 basis%bf = 0._dp
632 basis%dbf = 0._dp
633 basis%ddbf = 0._dp
634 DO l = 0, lmat
635 DO i = 1, basis%nbas(l)
636 al = basis%am(i, l)
637 DO k = 1, nr
638 rk = basis%grid%rad(k)
639 ear = exp(-al*basis%grid%rad(k)**2)
640 basis%bf(k, i, l) = rk**l*ear
641 basis%dbf(k, i, l) = (real(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear
642 basis%ddbf(k, i, l) = (real(l*(l - 1), dp)*rk**(l - 2) - &
643 2._dp*al*real(2*l + 1, dp)*rk**(l) + 4._dp*al*rk**(l + 2))*ear
644 END DO
645 END DO
646 END DO
647
648 CALL set_atom(atom, basis=basis)
649
650 ! optimization defaults
651 atom%optimization%damping = 0.2_dp
652 atom%optimization%eps_scf = 1.e-6_dp
653 atom%optimization%eps_diis = 100._dp
654 atom%optimization%max_iter = 50
655 atom%optimization%n_diis = 5
656
657 nelem = 0
658 ncore = 0
659 ncalc = 0
660 IF (monovalent) THEN
661 ncalc(0, 1) = 1
662 nelem(0, 1) = 1
663 ELSE IF (ASSOCIATED(gth_potential)) THEN
664 CALL get_potential(gth_potential, elec_conf=econf)
665 CALL set_pseudo_state(econf, z, ncalc, ncore, nelem)
666 ELSE IF (ASSOCIATED(sgp_potential)) THEN
667 CALL get_potential(sgp_potential, elec_conf=econf)
668 CALL set_pseudo_state(econf, z, ncalc, ncore, nelem)
669 ELSE
670 DO l = 0, min(lmat, ubound(ptable(z)%e_conv, 1))
671 ll = 2*(2*l + 1)
672 nn = ptable(z)%e_conv(l)
673 ii = 0
674 DO
675 ii = ii + 1
676 IF (nn <= ll) THEN
677 nelem(l, ii) = nn
678 EXIT
679 ELSE
680 nelem(l, ii) = ll
681 nn = nn - ll
682 END IF
683 END DO
684 END DO
685 ncalc = nelem - ncore
686 END IF
687
688 IF (qs_kind%ghost .OR. qs_kind%floating) THEN
689 nelem = 0
690 ncore = 0
691 ncalc = 0
692 END IF
693
694 ALLOCATE (atom%state)
695
696 atom%state%core = 0._dp
697 atom%state%core(0:lmat, 1:7) = real(ncore(0:lmat, 1:7), dp)
698 atom%state%occ = 0._dp
699 atom%state%occ(0:lmat, 1:7) = real(ncalc(0:lmat, 1:7), dp)
700 atom%state%occupation = 0._dp
701 atom%state%multiplicity = -1
702 DO l = 0, lmat
703 k = 0
704 DO i = 1, 7
705 IF (ncalc(l, i) > 0) THEN
706 k = k + 1
707 atom%state%occupation(l, k) = real(ncalc(l, i), dp)
708 END IF
709 END DO
710 END DO
711
712 atom%state%maxl_occ = get_maxl_occ(atom%state%occupation)
713 atom%state%maxn_occ = get_maxn_occ(atom%state%occupation)
714 atom%state%maxl_calc = atom%state%maxl_occ
715 atom%state%maxn_calc = atom%state%maxn_occ
716
717 ! calculate integrals
718 ! general integrals
719 CALL atom_int_setup(integrals, basis, potential=atom%potential, &
720 eri_coulomb=(atom%coulomb_integral_type == do_analytic), &
721 eri_exchange=(atom%exchange_integral_type == do_analytic))
722 ! potential
723 CALL atom_ppint_setup(integrals, basis, potential=atom%potential)
724 ! relativistic correction terms
725 NULLIFY (integrals%tzora, integrals%hdkh)
726 CALL atom_relint_setup(integrals, basis, atom%relativistic, zcore=real(atom%zcore, dp))
727 CALL set_atom(atom, integrals=integrals)
728
729 NULLIFY (orbitals)
730 mo = maxval(atom%state%maxn_calc)
731 mb = maxval(atom%basis%nbas)
732 CALL create_atom_orbs(orbitals, mb, mo)
733 CALL set_atom(atom, orbitals=orbitals)
734
735 CALL calculate_atom(atom, iw)
736
737 IF (do_basopt) THEN
738 CALL atom_fit_density(atom, ngto, 0, iw, results=results)
739 xx = results(1)
740 cc = results(2)
741 DO i = 1, ngto
742 density(i, 1) = xx*cc**i
743 density(i, 2) = results(2 + i)
744 END DO
745 ELSE
746 CALL atom_fit_density(atom, ngto, 0, iw, agto=density(:, 1), results=results)
747 density(1:ngto, 2) = results(1:ngto)
748 END IF
749
750 ! clean up
751 CALL atom_int_release(integrals)
752 CALL atom_ppint_release(integrals)
753 CALL atom_relint_release(integrals)
754 CALL release_atom_basis(basis)
755 CALL release_atom_potential(potential)
756 CALL release_atom_type(atom)
757
758 DEALLOCATE (potential, basis, integrals)
759
760 END SUBROUTINE calculate_atomic_density
761
762! **************************************************************************************************
763!> \brief ...
764!> \param atomic_kind ...
765!> \param qs_kind ...
766!> \param rel_control ...
767!> \param rtmat ...
768! **************************************************************************************************
769 SUBROUTINE calculate_atomic_relkin(atomic_kind, qs_kind, rel_control, rtmat)
770 TYPE(atomic_kind_type), INTENT(IN) :: atomic_kind
771 TYPE(qs_kind_type), INTENT(IN) :: qs_kind
772 TYPE(rel_control_type), POINTER :: rel_control
773 REAL(kind=dp), DIMENSION(:, :), POINTER :: rtmat
774
775 INTEGER :: i, ii, ipgf, j, k, k1, k2, l, ll, m, n, &
776 ngp, nj, nn, nr, ns, nset, nsgf, &
777 quadtype, relativistic, z
778 INTEGER, DIMENSION(0:lmat, 10) :: ncalc, ncore, nelem
779 INTEGER, DIMENSION(0:lmat, 100) :: set_index, shell_index
780 INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, nshell
781 INTEGER, DIMENSION(:, :), POINTER :: first_sgf, last_sgf, ls
782 REAL(kind=dp) :: al, alpha, ear, prefac, rk, zeff
783 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: omat
784 REAL(kind=dp), DIMENSION(:, :), POINTER :: zet
785 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: gcc
786 TYPE(all_potential_type), POINTER :: all_potential
787 TYPE(atom_basis_type), POINTER :: basis
788 TYPE(atom_integrals), POINTER :: integrals
789 TYPE(atom_potential_type), POINTER :: potential
790 TYPE(atom_type), POINTER :: atom
791 TYPE(grid_atom_type), POINTER :: grid
792 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
793
794 IF (rel_control%rel_method == rel_none) RETURN
795
796 NULLIFY (all_potential, orb_basis_set)
797 CALL get_qs_kind(qs_kind, basis_set=orb_basis_set, all_potential=all_potential)
798
799 cpassert(ASSOCIATED(orb_basis_set))
800
801 IF (ASSOCIATED(all_potential)) THEN
802 ! only all electron atoms will get the relativistic correction
803
804 CALL get_atomic_kind(atomic_kind, z=z)
805 CALL get_qs_kind(qs_kind, zeff=zeff)
806 NULLIFY (atom)
807 CALL create_atom_type(atom)
808 NULLIFY (atom%xc_section)
809 NULLIFY (atom%orbitals)
810 atom%z = z
811 alpha = sqrt(all_potential%alpha_core_charge)
812
813 ! set the method flag
814 SELECT CASE (rel_control%rel_method)
815 CASE DEFAULT
816 cpabort("")
817 CASE (rel_dkh)
818 SELECT CASE (rel_control%rel_DKH_order)
819 CASE DEFAULT
820 cpabort("")
821 CASE (0)
822 relativistic = do_dkh0_atom
823 CASE (1)
824 relativistic = do_dkh1_atom
825 CASE (2)
826 relativistic = do_dkh2_atom
827 CASE (3)
828 relativistic = do_dkh3_atom
829 END SELECT
830 CASE (rel_zora)
831 SELECT CASE (rel_control%rel_zora_type)
832 CASE DEFAULT
833 cpabort("")
834 CASE (rel_zora_full)
835 cpabort("")
836 CASE (rel_zora_mp)
837 relativistic = do_zoramp_atom
838 CASE (rel_sczora_mp)
839 relativistic = do_sczoramp_atom
840 END SELECT
841 END SELECT
842
843 CALL set_atom(atom, &
844 pp_calc=.false., &
845 method_type=do_rks_atom, &
846 relativistic=relativistic, &
847 coulomb_integral_type=do_numeric, &
848 exchange_integral_type=do_numeric)
849
850 ALLOCATE (potential, basis, integrals)
851
852 potential%ppot_type = no_pseudo
853 CALL set_atom(atom, zcore=z, potential=potential)
854
855 CALL get_gto_basis_set(orb_basis_set, &
856 nset=nset, nshell=nshell, npgf=npgf, lmin=lmin, lmax=lmax, l=ls, nsgf=nsgf, zet=zet, gcc=gcc, &
857 first_sgf=first_sgf, last_sgf=last_sgf)
858
859 NULLIFY (grid)
860 ngp = 400
861 quadtype = do_gapw_log
862 CALL allocate_grid_atom(grid)
863 CALL create_grid_atom(grid, ngp, 1, 1, 0, quadtype)
864 grid%nr = ngp
865 basis%grid => grid
866
867 NULLIFY (basis%am, basis%cm, basis%as, basis%ns, basis%bf, basis%dbf, basis%ddbf)
868 basis%basis_type = cgto_basis
869 basis%eps_eig = 1.e-12_dp
870
871 ! fill in the basis data structures
872 set_index = 0
873 shell_index = 0
874 basis%nprim = 0
875 basis%nbas = 0
876 DO i = 1, nset
877 DO j = lmin(i), min(lmax(i), lmat)
878 basis%nprim(j) = basis%nprim(j) + npgf(i)
879 END DO
880 DO j = 1, nshell(i)
881 l = ls(j, i)
882 IF (l <= lmat) THEN
883 basis%nbas(l) = basis%nbas(l) + 1
884 k = basis%nbas(l)
885 cpassert(k <= 100)
886 set_index(l, k) = i
887 shell_index(l, k) = j
888 END IF
889 END DO
890 END DO
891
892 nj = maxval(basis%nprim)
893 ns = maxval(basis%nbas)
894 ALLOCATE (basis%am(nj, 0:lmat))
895 basis%am = 0._dp
896 ALLOCATE (basis%cm(nj, ns, 0:lmat))
897 basis%cm = 0._dp
898 DO j = 0, lmat
899 nj = 0
900 ns = 0
901 DO i = 1, nset
902 IF (j >= lmin(i) .AND. j <= lmax(i)) THEN
903 DO ipgf = 1, npgf(i)
904 basis%am(nj + ipgf, j) = zet(ipgf, i)
905 END DO
906 DO ii = 1, nshell(i)
907 IF (ls(ii, i) == j) THEN
908 ns = ns + 1
909 DO ipgf = 1, npgf(i)
910 basis%cm(nj + ipgf, ns, j) = gcc(ipgf, ii, i)
911 END DO
912 END IF
913 END DO
914 nj = nj + npgf(i)
915 END IF
916 END DO
917 END DO
918
919 ! Normalization as used in the atomic code
920 ! We have to undo the Quickstep normalization
921 DO j = 0, lmat
922 prefac = 2.0_dp*sqrt(pi/dfac(2*j + 1))
923 DO ipgf = 1, basis%nprim(j)
924 DO ii = 1, basis%nbas(j)
925 basis%cm(ipgf, ii, j) = prefac*basis%cm(ipgf, ii, j)
926 END DO
927 END DO
928 END DO
929
930 ! initialize basis function on a radial grid
931 nr = basis%grid%nr
932 m = maxval(basis%nbas)
933 ALLOCATE (basis%bf(nr, m, 0:lmat))
934 ALLOCATE (basis%dbf(nr, m, 0:lmat))
935 ALLOCATE (basis%ddbf(nr, m, 0:lmat))
936
937 basis%bf = 0._dp
938 basis%dbf = 0._dp
939 basis%ddbf = 0._dp
940 DO l = 0, lmat
941 DO i = 1, basis%nprim(l)
942 al = basis%am(i, l)
943 DO k = 1, nr
944 rk = basis%grid%rad(k)
945 ear = exp(-al*basis%grid%rad(k)**2)
946 DO j = 1, basis%nbas(l)
947 basis%bf(k, j, l) = basis%bf(k, j, l) + rk**l*ear*basis%cm(i, j, l)
948 basis%dbf(k, j, l) = basis%dbf(k, j, l) &
949 + (real(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear*basis%cm(i, j, l)
950 basis%ddbf(k, j, l) = basis%ddbf(k, j, l) + &
951 (real(l*(l - 1), dp)*rk**(l - 2) - 2._dp*al*real(2*l + 1, dp)* &
952 rk**(l) + 4._dp*al*rk**(l + 2))*ear*basis%cm(i, j, l)
953 END DO
954 END DO
955 END DO
956 END DO
957
958 CALL set_atom(atom, basis=basis)
959
960 ! optimization defaults
961 atom%optimization%damping = 0.2_dp
962 atom%optimization%eps_scf = 1.e-6_dp
963 atom%optimization%eps_diis = 100._dp
964 atom%optimization%max_iter = 50
965 atom%optimization%n_diis = 5
966
967 ! electronic state
968 nelem = 0
969 ncore = 0
970 ncalc = 0
971 DO l = 0, min(lmat, ubound(ptable(z)%e_conv, 1))
972 ll = 2*(2*l + 1)
973 nn = ptable(z)%e_conv(l)
974 ii = 0
975 DO
976 ii = ii + 1
977 IF (nn <= ll) THEN
978 nelem(l, ii) = nn
979 EXIT
980 ELSE
981 nelem(l, ii) = ll
982 nn = nn - ll
983 END IF
984 END DO
985 END DO
986 ncalc = nelem - ncore
987
988 IF (qs_kind%ghost .OR. qs_kind%floating) THEN
989 nelem = 0
990 ncore = 0
991 ncalc = 0
992 END IF
993
994 ALLOCATE (atom%state)
995
996 atom%state%core = 0._dp
997 atom%state%core(0:lmat, 1:7) = real(ncore(0:lmat, 1:7), dp)
998 atom%state%occ = 0._dp
999 atom%state%occ(0:lmat, 1:7) = real(ncalc(0:lmat, 1:7), dp)
1000 atom%state%occupation = 0._dp
1001 atom%state%multiplicity = -1
1002 DO l = 0, lmat
1003 k = 0
1004 DO i = 1, 7
1005 IF (ncalc(l, i) > 0) THEN
1006 k = k + 1
1007 atom%state%occupation(l, k) = real(ncalc(l, i), dp)
1008 END IF
1009 END DO
1010 END DO
1011
1012 atom%state%maxl_occ = get_maxl_occ(atom%state%occupation)
1013 atom%state%maxn_occ = get_maxn_occ(atom%state%occupation)
1014 atom%state%maxl_calc = atom%state%maxl_occ
1015 atom%state%maxn_calc = atom%state%maxn_occ
1016
1017 ! calculate integrals
1018 ! general integrals
1019 CALL atom_int_setup(integrals, basis)
1020 ! potential
1021 CALL atom_ppint_setup(integrals, basis, potential=atom%potential)
1022 ! relativistic correction terms
1023 NULLIFY (integrals%tzora, integrals%hdkh)
1024 CALL atom_relint_setup(integrals, basis, atom%relativistic, zcore=real(atom%zcore, dp), &
1025 alpha=alpha)
1026 CALL set_atom(atom, integrals=integrals)
1027
1028 ! for DKH we need erfc integrals to correct non-relativistic
1029 integrals%core = 0.0_dp
1030 DO l = 0, lmat
1031 n = integrals%n(l)
1032 m = basis%nprim(l)
1033 ALLOCATE (omat(m, m))
1034
1035 CALL sg_erfc(omat(1:m, 1:m), l, alpha, basis%am(1:m, l), basis%am(1:m, l))
1036 integrals%core(1:n, 1:n, l) = matmul(transpose(basis%cm(1:m, 1:n, l)), &
1037 matmul(omat(1:m, 1:m), basis%cm(1:m, 1:n, l)))
1038
1039 DEALLOCATE (omat)
1040 END DO
1041
1042 ! recover relativistic kinetic matrix in CP2K/GPW order and normalization
1043 IF (ASSOCIATED(rtmat)) THEN
1044 DEALLOCATE (rtmat)
1045 END IF
1046 ALLOCATE (rtmat(nsgf, nsgf))
1047 rtmat = 0._dp
1048 DO l = 0, lmat
1049 ll = 2*l
1050 DO k1 = 1, basis%nbas(l)
1051 DO k2 = 1, basis%nbas(l)
1052 i = first_sgf(shell_index(l, k1), set_index(l, k1))
1053 j = first_sgf(shell_index(l, k2), set_index(l, k2))
1054 SELECT CASE (atom%relativistic)
1055 CASE DEFAULT
1056 cpabort("")
1057 CASE (do_zoramp_atom, do_sczoramp_atom)
1058 DO m = 0, ll
1059 rtmat(i + m, j + m) = integrals%tzora(k1, k2, l)
1060 END DO
1061 CASE (do_dkh0_atom, do_dkh1_atom, do_dkh2_atom, do_dkh3_atom)
1062 DO m = 0, ll
1063 rtmat(i + m, j + m) = integrals%hdkh(k1, k2, l) - integrals%kin(k1, k2, l) + &
1064 atom%zcore*integrals%core(k1, k2, l)
1065 END DO
1066 END SELECT
1067 END DO
1068 END DO
1069 END DO
1070 DO k1 = 1, nsgf
1071 DO k2 = k1, nsgf
1072 rtmat(k1, k2) = 0.5_dp*(rtmat(k1, k2) + rtmat(k2, k1))
1073 rtmat(k2, k1) = rtmat(k1, k2)
1074 END DO
1075 END DO
1076
1077 ! clean up
1078 CALL atom_int_release(integrals)
1079 CALL atom_ppint_release(integrals)
1080 CALL atom_relint_release(integrals)
1081 CALL release_atom_basis(basis)
1082 CALL release_atom_potential(potential)
1083 CALL release_atom_type(atom)
1084
1085 DEALLOCATE (potential, basis, integrals)
1086
1087 ELSE
1088
1089 IF (ASSOCIATED(rtmat)) THEN
1090 DEALLOCATE (rtmat)
1091 END IF
1092 NULLIFY (rtmat)
1093
1094 END IF
1095
1096 END SUBROUTINE calculate_atomic_relkin
1097
1098! **************************************************************************************************
1099!> \brief ...
1100!> \param gth_potential ...
1101!> \param gth_atompot ...
1102! **************************************************************************************************
1103 SUBROUTINE gth_potential_conversion(gth_potential, gth_atompot)
1104 TYPE(gth_potential_type), POINTER :: gth_potential
1105 TYPE(atom_gthpot_type) :: gth_atompot
1106
1107 INTEGER :: i, j, l, lm, n, ne, nexp_lpot, nexp_lsd, &
1108 nexp_nlcc
1109 INTEGER, DIMENSION(:), POINTER :: nct_lpot, nct_lsd, nct_nlcc, nppnl, &
1110 ppeconf
1111 LOGICAL :: lpot_present, lsd_present, nlcc_present, &
1112 soc_present
1113 REAL(kind=dp) :: ac, zeff
1114 REAL(kind=dp), DIMENSION(:), POINTER :: alpha_lpot, alpha_lsd, alpha_nlcc, ap, ce
1115 REAL(kind=dp), DIMENSION(:, :), POINTER :: cval_lpot, cval_lsd, cval_nlcc
1116 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: hp, kp
1117
1118 CALL get_potential(gth_potential, &
1119 zeff=zeff, &
1120 elec_conf=ppeconf, &
1121 alpha_core_charge=ac, &
1122 nexp_ppl=ne, &
1123 cexp_ppl=ce, &
1124 lppnl=lm, &
1125 nprj_ppnl=nppnl, &
1126 alpha_ppnl=ap, &
1127 kprj_ppnl=kp, &
1128 hprj_ppnl=hp)
1129
1130 gth_atompot%zion = zeff
1131 gth_atompot%rc = sqrt(0.5_dp/ac)
1132 gth_atompot%ncl = ne
1133 gth_atompot%cl(:) = 0._dp
1134 IF (ac > 0._dp) THEN
1135 DO i = 1, ne
1136 gth_atompot%cl(i) = ce(i)/(2._dp*ac)**(i - 1)
1137 END DO
1138 END IF
1139 !extended type
1140 gth_atompot%lpotextended = .false.
1141 gth_atompot%lsdpot = .false.
1142 gth_atompot%nlcc = .false.
1143 gth_atompot%nexp_lpot = 0
1144 gth_atompot%nexp_lsd = 0
1145 gth_atompot%nexp_nlcc = 0
1146 CALL get_potential(gth_potential, &
1147 lpot_present=lpot_present, &
1148 lsd_present=lsd_present, &
1149 nlcc_present=nlcc_present)
1150 IF (lpot_present) THEN
1151 CALL get_potential(gth_potential, &
1152 nexp_lpot=nexp_lpot, &
1153 alpha_lpot=alpha_lpot, &
1154 nct_lpot=nct_lpot, &
1155 cval_lpot=cval_lpot)
1156 gth_atompot%lpotextended = .true.
1157 gth_atompot%nexp_lpot = nexp_lpot
1158 gth_atompot%alpha_lpot(1:nexp_lpot) = sqrt(0.5_dp/alpha_lpot(1:nexp_lpot))
1159 gth_atompot%nct_lpot(1:nexp_lpot) = nct_lpot(1:nexp_lpot)
1160 DO j = 1, nexp_lpot
1161 ac = alpha_lpot(j)
1162 DO i = 1, 4
1163 gth_atompot%cval_lpot(i, j) = cval_lpot(i, j)/(2._dp*ac)**(i - 1)
1164 END DO
1165 END DO
1166 END IF
1167 IF (lsd_present) THEN
1168 CALL get_potential(gth_potential, &
1169 nexp_lsd=nexp_lsd, &
1170 alpha_lsd=alpha_lsd, &
1171 nct_lsd=nct_lsd, &
1172 cval_lsd=cval_lsd)
1173 gth_atompot%lsdpot = .true.
1174 gth_atompot%nexp_lsd = nexp_lsd
1175 gth_atompot%alpha_lsd(1:nexp_lsd) = sqrt(0.5_dp/alpha_lsd(1:nexp_lsd))
1176 gth_atompot%nct_lsd(1:nexp_lsd) = nct_lsd(1:nexp_lsd)
1177 DO j = 1, nexp_lpot
1178 ac = alpha_lsd(j)
1179 DO i = 1, 4
1180 gth_atompot%cval_lsd(i, j) = cval_lsd(i, j)/(2._dp*ac)**(i - 1)
1181 END DO
1182 END DO
1183 END IF
1184
1185 ! nonlocal part
1186 gth_atompot%nl(:) = 0
1187 gth_atompot%rcnl(:) = 0._dp
1188 gth_atompot%hnl(:, :, :) = 0._dp
1189 DO l = 0, lm
1190 n = nppnl(l)
1191 gth_atompot%nl(l) = n
1192 gth_atompot%rcnl(l) = sqrt(0.5_dp/ap(l))
1193 gth_atompot%hnl(1:n, 1:n, l) = hp(1:n, 1:n, l)
1194 END DO
1195
1196 ! SOC
1197 CALL get_potential(gth_potential, soc_present=soc_present)
1198 gth_atompot%soc = soc_present
1199 gth_atompot%knl = 0.0_dp
1200 IF (soc_present) THEN
1201 DO l = 1, lm
1202 n = nppnl(l)
1203 gth_atompot%knl(1:n, 1:n, l) = kp(1:n, 1:n, l)
1204 END DO
1205 END IF
1206
1207 IF (nlcc_present) THEN
1208 CALL get_potential(gth_potential, &
1209 nexp_nlcc=nexp_nlcc, &
1210 alpha_nlcc=alpha_nlcc, &
1211 nct_nlcc=nct_nlcc, &
1212 cval_nlcc=cval_nlcc)
1213 gth_atompot%nlcc = .true.
1214 gth_atompot%nexp_nlcc = nexp_nlcc
1215 gth_atompot%alpha_nlcc(1:nexp_nlcc) = alpha_nlcc(1:nexp_nlcc)
1216 gth_atompot%nct_nlcc(1:nexp_nlcc) = nct_nlcc(1:nexp_nlcc)
1217 gth_atompot%cval_nlcc(1:4, 1:nexp_nlcc) = cval_nlcc(1:4, 1:nexp_nlcc)
1218 END IF
1219
1220 END SUBROUTINE gth_potential_conversion
1221
1222! **************************************************************************************************
1223!> \brief ...
1224!> \param sgp_potential ...
1225!> \param sgp_atompot ...
1226! **************************************************************************************************
1227 SUBROUTINE sgp_potential_conversion(sgp_potential, sgp_atompot)
1228 TYPE(sgp_potential_type), POINTER :: sgp_potential
1229 TYPE(atom_sgppot_type) :: sgp_atompot
1230
1231 INTEGER :: lm, n
1232 INTEGER, DIMENSION(:), POINTER :: ppeconf
1233 LOGICAL :: nlcc_present
1234 REAL(kind=dp) :: ac, zeff
1235 REAL(kind=dp), DIMENSION(:), POINTER :: ap, ce
1236 REAL(kind=dp), DIMENSION(:, :), POINTER :: hhp
1237 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: ccp
1238
1239 CALL get_potential(sgp_potential, &
1240 name=sgp_atompot%pname, &
1241 zeff=zeff, &
1242 elec_conf=ppeconf, &
1243 alpha_core_charge=ac)
1244 sgp_atompot%zion = zeff
1245 sgp_atompot%ac_local = ac
1246 sgp_atompot%econf(0:3) = ppeconf(0:3)
1247 CALL get_potential(sgp_potential, lmax=lm, &
1248 is_nonlocal=sgp_atompot%is_nonlocal, &
1249 n_nonlocal=n, a_nonlocal=ap, h_nonlocal=hhp, c_nonlocal=ccp)
1250 ! nonlocal
1251 sgp_atompot%has_nonlocal = any(sgp_atompot%is_nonlocal)
1252 sgp_atompot%lmax = lm
1253 IF (sgp_atompot%has_nonlocal) THEN
1254 cpassert(n <= SIZE(sgp_atompot%a_nonlocal))
1255 sgp_atompot%n_nonlocal = n
1256 sgp_atompot%a_nonlocal(1:n) = ap(1:n)
1257 sgp_atompot%h_nonlocal(1:n, 0:lm) = hhp(1:n, 0:lm)
1258 sgp_atompot%c_nonlocal(1:n, 1:n, 0:lm) = ccp(1:n, 1:n, 0:lm)
1259 END IF
1260 ! local
1261 CALL get_potential(sgp_potential, n_local=n, a_local=ap, c_local=ce)
1262 cpassert(n <= SIZE(sgp_atompot%a_local))
1263 sgp_atompot%n_local = n
1264 sgp_atompot%a_local(1:n) = ap(1:n)
1265 sgp_atompot%c_local(1:n) = ce(1:n)
1266 ! NLCC
1267 CALL get_potential(sgp_potential, has_nlcc=nlcc_present, &
1268 n_nlcc=n, a_nlcc=ap, c_nlcc=ce)
1269 IF (nlcc_present) THEN
1270 sgp_atompot%has_nlcc = .true.
1271 cpassert(n <= SIZE(sgp_atompot%a_nlcc))
1272 sgp_atompot%n_nlcc = n
1273 sgp_atompot%a_nlcc(1:n) = ap(1:n)
1274 sgp_atompot%c_nlcc(1:n) = ce(1:n)
1275 ELSE
1276 sgp_atompot%has_nlcc = .false.
1277 END IF
1278
1279 END SUBROUTINE sgp_potential_conversion
1280
1281! **************************************************************************************************
1282!> \brief ...
1283!> \param sgp_potential ...
1284!> \param ecp_atompot ...
1285! **************************************************************************************************
1286 SUBROUTINE ecp_potential_conversion(sgp_potential, ecp_atompot)
1287 TYPE(sgp_potential_type), POINTER :: sgp_potential
1288 TYPE(atom_ecppot_type) :: ecp_atompot
1289
1290 INTEGER, DIMENSION(:), POINTER :: ppeconf
1291 LOGICAL :: ecp_local, ecp_semi_local
1292 REAL(kind=dp) :: zeff
1293
1294 CALL get_potential(sgp_potential, ecp_local=ecp_local, ecp_semi_local=ecp_semi_local)
1295 cpassert(ecp_semi_local .AND. ecp_local)
1296 CALL get_potential(sgp_potential, &
1297 name=ecp_atompot%pname, &
1298 zeff=zeff, &
1299 elec_conf=ppeconf)
1300 ecp_atompot%zion = zeff
1301 ecp_atompot%econf(0:3) = ppeconf(0:3)
1302 CALL get_potential(sgp_potential, sl_lmax=ecp_atompot%lmax)
1303 ! local
1304 CALL get_potential(sgp_potential, nloc=ecp_atompot%nloc, nrloc=ecp_atompot%nrloc, &
1305 aloc=ecp_atompot%aloc, bloc=ecp_atompot%bloc)
1306 ! nonlocal
1307 CALL get_potential(sgp_potential, npot=ecp_atompot%npot, nrpot=ecp_atompot%nrpot, &
1308 apot=ecp_atompot%apot, bpot=ecp_atompot%bpot)
1309
1310 END SUBROUTINE ecp_potential_conversion
1311! **************************************************************************************************
1312
1313END MODULE atom_kind_orbitals
external_potential_types::get_potential
Definition external_potential_types.F:278
ai_onecenter
Definition ai_onecenter.F:28
ai_onecenter::sg_erfc
subroutine, public sg_erfc(umat, l, a, pa, pb)
...
Definition ai_onecenter.F:590
atom_electronic_structure
Definition atom_electronic_structure.F:8
atom_electronic_structure::calculate_atom
subroutine, public calculate_atom(atom, iw, noguess, converged)
General routine to perform electronic structure atomic calculations.
Definition atom_electronic_structure.F:60
atom_fit
routines that fit parameters for /from atomic calculations
Definition atom_fit.F:11
atom_fit::atom_fit_density
subroutine, public atom_fit_density(atom, num_gto, norder, iunit, agto, powell_section, results)
Fit the atomic electron density using a geometrical Gaussian basis set.
Definition atom_fit.F:78
atom_kind_orbitals
calculate the orbitals for a given atomic kind type
Definition atom_kind_orbitals.F:11
atom_kind_orbitals::calculate_atomic_relkin
subroutine, public calculate_atomic_relkin(atomic_kind, qs_kind, rel_control, rtmat)
...
Definition atom_kind_orbitals.F:770
atom_kind_orbitals::calculate_atomic_orbitals
subroutine, public calculate_atomic_orbitals(atomic_kind, qs_kind, agrid, iunit, pmat, fmat, density, wavefunction, wfninfo, confine, xc_section, nocc)
...
Definition atom_kind_orbitals.F:89
atom_kind_orbitals::calculate_atomic_density
subroutine, public calculate_atomic_density(density, atomic_kind, qs_kind, ngto, iunit, optbasis, allelectron, confine)
...
Definition atom_kind_orbitals.F:481
atom_kind_orbitals::gth_potential_conversion
subroutine, public gth_potential_conversion(gth_potential, gth_atompot)
...
Definition atom_kind_orbitals.F:1104
atom_operators
Calculate the atomic operator matrices.
Definition atom_operators.F:15
atom_operators::atom_ppint_release
subroutine, public atom_ppint_release(integrals)
Release memory allocated for atomic integrals (core electrons).
Definition atom_operators.F:1160
atom_operators::atom_int_setup
subroutine, public atom_int_setup(integrals, basis, potential, eri_coulomb, eri_exchange, all_nu)
Set up atomic integrals.
Definition atom_operators.F:63
atom_operators::atom_relint_setup
subroutine, public atom_relint_setup(integrals, basis, reltyp, zcore, alpha)
...
Definition atom_operators.F:817
atom_operators::atom_relint_release
subroutine, public atom_relint_release(integrals)
Release memory allocated for atomic integrals (relativistic effects).
Definition atom_operators.F:1181
atom_operators::atom_ppint_setup
subroutine, public atom_ppint_setup(integrals, basis, potential)
...
Definition atom_operators.F:336
atom_operators::atom_int_release
subroutine, public atom_int_release(integrals)
Release memory allocated for atomic integrals (valence electrons).
Definition atom_operators.F:1123
atom_set_basis
Definition atom_set_basis.F:9
atom_set_basis::set_kind_basis_atomic
subroutine, public set_kind_basis_atomic(basis, orb_basis_set, has_pp, agrid, cp2k_norm)
...
Definition atom_set_basis.F:58
atom_types
Define the atom type and its sub types.
Definition atom_types.F:15
atom_types::create_atom_type
subroutine, public create_atom_type(atom)
...
Definition atom_types.F:958
atom_types::cgto_basis
integer, parameter, public cgto_basis
Definition atom_types.F:69
atom_types::gto_basis
integer, parameter, public gto_basis
Definition atom_types.F:69
atom_types::release_atom_type
subroutine, public release_atom_type(atom)
...
Definition atom_types.F:982
atom_types::release_atom_potential
subroutine, public release_atom_potential(potential)
...
Definition atom_types.F:2535
atom_types::lmat
integer, parameter, public lmat
Definition atom_types.F:67
atom_types::set_atom
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)
...
Definition atom_types.F:1023
atom_types::release_atom_basis
subroutine, public release_atom_basis(basis)
...
Definition atom_types.F:924
atom_types::create_atom_orbs
subroutine, public create_atom_orbs(orbs, mbas, mo)
...
Definition atom_types.F:1069
atom_types::clementi_geobas
subroutine, public clementi_geobas(zval, cval, aval, ngto, ival)
...
Definition atom_types.F:1442
atom_utils
Some basic routines for atomic calculations.
Definition atom_utils.F:15
atom_utils::get_maxn_occ
pure integer function, dimension(0:lmat), public get_maxn_occ(occupation)
Return the maximum principal quantum number of occupied orbitals.
Definition atom_utils.F:301
atom_utils::atom_density
subroutine, public atom_density(density, pmat, basis, maxl, typ, rr)
Map the electron density on an atomic radial grid.
Definition atom_utils.F:366
atom_utils::get_maxl_occ
pure integer function, public get_maxl_occ(occupation)
Return the maximum orbital quantum number of occupied orbitals.
Definition atom_utils.F:281
atom
Definition atom.F:9
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
basis_set_types
Definition basis_set_types.F:15
basis_set_types::get_gto_basis_set
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum)
...
Definition basis_set_types.F:631
external_potential_types
Definition of the atomic potential types.
Definition external_potential_types.F:14
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::rel_zora_full
integer, parameter, public rel_zora_full
Definition input_constants.F:612
input_constants::do_rks_atom
integer, parameter, public do_rks_atom
Definition input_constants.F:159
input_constants::do_analytic
integer, parameter, public do_analytic
Definition input_constants.F:149
input_constants::sgp_pseudo
integer, parameter, public sgp_pseudo
Definition input_constants.F:174
input_constants::do_dkh3_atom
integer, parameter, public do_dkh3_atom
Definition input_constants.F:165
input_constants::gth_pseudo
integer, parameter, public gth_pseudo
Definition input_constants.F:174
input_constants::ecp_pseudo
integer, parameter, public ecp_pseudo
Definition input_constants.F:174
input_constants::do_nonrel_atom
integer, parameter, public do_nonrel_atom
Definition input_constants.F:165
input_constants::do_dkh0_atom
integer, parameter, public do_dkh0_atom
Definition input_constants.F:165
input_constants::rel_zora_mp
integer, parameter, public rel_zora_mp
Definition input_constants.F:612
input_constants::rel_zora
integer, parameter, public rel_zora
Definition input_constants.F:609
input_constants::poly_conf
integer, parameter, public poly_conf
Definition input_constants.F:180
input_constants::do_dkh2_atom
integer, parameter, public do_dkh2_atom
Definition input_constants.F:165
input_constants::no_pseudo
integer, parameter, public no_pseudo
Definition input_constants.F:174
input_constants::do_uks_atom
integer, parameter, public do_uks_atom
Definition input_constants.F:159
input_constants::barrier_conf
integer, parameter, public barrier_conf
Definition input_constants.F:180
input_constants::rel_dkh
integer, parameter, public rel_dkh
Definition input_constants.F:609
input_constants::do_numeric
integer, parameter, public do_numeric
Definition input_constants.F:149
input_constants::do_zoramp_atom
integer, parameter, public do_zoramp_atom
Definition input_constants.F:165
input_constants::do_gapw_log
integer, parameter, public do_gapw_log
Definition input_constants.F:226
input_constants::do_dkh1_atom
integer, parameter, public do_dkh1_atom
Definition input_constants.F:165
input_constants::rel_none
integer, parameter, public rel_none
Definition input_constants.F:609
input_constants::rel_sczora_mp
integer, parameter, public rel_sczora_mp
Definition input_constants.F:612
input_constants::do_sczoramp_atom
integer, parameter, public do_sczoramp_atom
Definition input_constants.F:165
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
mathconstants::dfac
real(kind=dp), dimension(-1:2 *maxfac+1), parameter, public dfac
Definition mathconstants.F:61
periodic_table
Periodic Table related data definitions.
Definition periodic_table.F:28
periodic_table::ptable
type(atom), dimension(0:nelem), public ptable
Definition periodic_table.F:65
periodic_table::nelem
integer, parameter, public nelem
Definition periodic_table.F:63
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::bohr
real(kind=dp), parameter, public bohr
Definition physcon.F:147
qs_grid_atom
Definition qs_grid_atom.F:8
qs_grid_atom::allocate_grid_atom
subroutine, public allocate_grid_atom(grid_atom)
Initialize components of the grid_atom_type structure.
Definition qs_grid_atom.F:73
qs_grid_atom::create_grid_atom
subroutine, public create_grid_atom(grid_atom, nr, na, llmax, ll, quadrature)
...
Definition qs_grid_atom.F:183
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::set_pseudo_state
subroutine, public set_pseudo_state(econf, z, ncalc, ncore, nelem)
...
Definition qs_kind_types.F:3547
qs_kind_types::has_nlcc
logical function, public has_nlcc(qs_kind_set)
finds if a given qs run needs to use nlcc
Definition qs_kind_types.F:3973
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:468
qs_kind_types::init_atom_electronic_state
subroutine, public init_atom_electronic_state(atomic_kind, qs_kind, ncalc, ncore, nelem, edelta)
...
Definition qs_kind_types.F:3410
rel_control_types
parameters that control a relativistic calculation
Definition rel_control_types.F:15
atom_types::atom_basis_type
Provides all information about a basis set.
Definition atom_types.F:78
atom_types::atom_ecppot_type
Definition atom_types.F:133
atom_types::atom_gthpot_type
Provides all information about a pseudopotential.
Definition atom_types.F:98
atom_types::atom_integrals
Definition atom_types.F:216
atom_types::atom_orbitals
Holds atomic orbitals and energies.
Definition atom_types.F:237
atom_types::atom_potential_type
Definition atom_types.F:171
atom_types::atom_sgppot_type
Definition atom_types.F:149
atom_types::atom_type
Provides all information about an atomic kind.
Definition atom_types.F:293
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
external_potential_types::all_potential_type
Definition external_potential_types.F:70
external_potential_types::gth_potential_type
Definition external_potential_types.F:116
external_potential_types::sgp_potential_type
Definition external_potential_types.F:174
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
qs_grid_atom::grid_atom_type
Definition qs_grid_atom.F:44
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:177
rel_control_types::rel_control_type
contains the parameters needed by a relativistic calculation
Definition rel_control_types.F:55