(git:374b731)
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-2024 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 allelectron ...
476!> \param confine ...
477! **************************************************************************************************
478 SUBROUTINE calculate_atomic_density(density, atomic_kind, qs_kind, ngto, iunit, allelectron, confine)
479 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: density
480 TYPE(atomic_kind_type), POINTER :: atomic_kind
481 TYPE(qs_kind_type), POINTER :: qs_kind
482 INTEGER, INTENT(IN) :: ngto
483 INTEGER, INTENT(IN), OPTIONAL :: iunit
484 LOGICAL, INTENT(IN), OPTIONAL :: allelectron, confine
485
486 INTEGER, PARAMETER :: num_gto = 40
487
488 INTEGER :: i, ii, k, l, ll, m, mb, mo, ngp, nn, nr, &
489 quadtype, relativistic, z
490 INTEGER, DIMENSION(0:lmat) :: starti
491 INTEGER, DIMENSION(0:lmat, 10) :: ncalc, ncore, nelem
492 INTEGER, DIMENSION(:), POINTER :: econf
493 LOGICAL :: ecp_semi_local
494 REAL(kind=dp) :: al, aval, cc, cval, ear, rk, xx, zeff
495 REAL(kind=dp), DIMENSION(num_gto+2) :: results
496 TYPE(all_potential_type), POINTER :: all_potential
497 TYPE(atom_basis_type), POINTER :: basis
498 TYPE(atom_integrals), POINTER :: integrals
499 TYPE(atom_orbitals), POINTER :: orbitals
500 TYPE(atom_potential_type), POINTER :: potential
501 TYPE(atom_type), POINTER :: atom
502 TYPE(grid_atom_type), POINTER :: grid
503 TYPE(gth_potential_type), POINTER :: gth_potential
504 TYPE(sgp_potential_type), POINTER :: sgp_potential
505
506 NULLIFY (atom)
507 CALL create_atom_type(atom)
508
509 CALL get_atomic_kind(atomic_kind, z=z)
510 NULLIFY (all_potential, gth_potential)
511 CALL get_qs_kind(qs_kind, zeff=zeff, &
512 all_potential=all_potential, &
513 gth_potential=gth_potential, &
514 sgp_potential=sgp_potential)
515
516 IF (PRESENT(allelectron)) THEN
517 IF (allelectron) THEN
518 NULLIFY (gth_potential)
519 zeff = z
520 END IF
521 END IF
522
523 cpassert(ngto <= num_gto)
524
525 IF (ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)) THEN
526 ! PP calculation are non-relativistic
527 relativistic = do_nonrel_atom
528 ELSE
529 ! AE calculations use DKH2
530 relativistic = do_dkh2_atom
531 END IF
532
533 atom%z = z
534 CALL set_atom(atom, &
535 pp_calc=(ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)), &
536 method_type=do_rks_atom, &
537 relativistic=relativistic, &
538 coulomb_integral_type=do_numeric, &
539 exchange_integral_type=do_numeric)
540
541 ALLOCATE (potential, basis, integrals)
542
543 IF (PRESENT(confine)) THEN
544 potential%confinement = confine
545 ELSE
546 IF (ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)) THEN
547 potential%confinement = .true.
548 ELSE
549 potential%confinement = .false.
550 END IF
551 END IF
552 potential%conf_type = barrier_conf
553 potential%acon = 200._dp
554 potential%rcon = 4.0_dp
555 potential%scon = 8.0_dp
556
557 IF (ASSOCIATED(gth_potential)) THEN
558 potential%ppot_type = gth_pseudo
559 CALL get_potential(gth_potential, zeff=zeff)
560 CALL gth_potential_conversion(gth_potential, potential%gth_pot)
561 CALL set_atom(atom, zcore=nint(zeff), potential=potential)
562 ELSE IF (ASSOCIATED(sgp_potential)) THEN
563 CALL get_potential(sgp_potential, ecp_semi_local=ecp_semi_local)
564 IF (ecp_semi_local) THEN
565 potential%ppot_type = ecp_pseudo
566 CALL ecp_potential_conversion(sgp_potential, potential%ecp_pot)
567 potential%ecp_pot%symbol = ptable(z)%symbol
568 ELSE
569 potential%ppot_type = sgp_pseudo
570 CALL sgp_potential_conversion(sgp_potential, potential%sgp_pot)
571 potential%sgp_pot%symbol = ptable(z)%symbol
572 END IF
573 CALL get_potential(sgp_potential, zeff=zeff)
574 CALL set_atom(atom, zcore=nint(zeff), potential=potential)
575 ELSE
576 potential%ppot_type = no_pseudo
577 CALL set_atom(atom, zcore=z, potential=potential)
578 END IF
579
580 ! atomic grid
581 NULLIFY (grid)
582 ngp = 400
583 quadtype = do_gapw_log
584 CALL allocate_grid_atom(grid)
585 CALL create_grid_atom(grid, ngp, 1, 1, 0, quadtype)
586 grid%nr = ngp
587 basis%grid => grid
588
589 NULLIFY (basis%am, basis%cm, basis%as, basis%ns, basis%bf, basis%dbf, basis%ddbf)
590
591 ! fill in the basis data structures
592 basis%eps_eig = 1.e-12_dp
593 basis%basis_type = gto_basis
594 CALL clementi_geobas(z, cval, aval, basis%nbas, starti)
595 basis%nprim = basis%nbas
596 m = maxval(basis%nbas)
597 ALLOCATE (basis%am(m, 0:lmat))
598 basis%am = 0._dp
599 DO l = 0, lmat
600 DO i = 1, basis%nbas(l)
601 ll = i - 1 + starti(l)
602 basis%am(i, l) = aval*cval**(ll)
603 END DO
604 END DO
605
606 basis%geometrical = .true.
607 basis%aval = aval
608 basis%cval = cval
609 basis%start = starti
610
611 ! initialize basis function on a radial grid
612 nr = basis%grid%nr
613 m = maxval(basis%nbas)
614 ALLOCATE (basis%bf(nr, m, 0:lmat))
615 ALLOCATE (basis%dbf(nr, m, 0:lmat))
616 ALLOCATE (basis%ddbf(nr, m, 0:lmat))
617 basis%bf = 0._dp
618 basis%dbf = 0._dp
619 basis%ddbf = 0._dp
620 DO l = 0, lmat
621 DO i = 1, basis%nbas(l)
622 al = basis%am(i, l)
623 DO k = 1, nr
624 rk = basis%grid%rad(k)
625 ear = exp(-al*basis%grid%rad(k)**2)
626 basis%bf(k, i, l) = rk**l*ear
627 basis%dbf(k, i, l) = (real(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear
628 basis%ddbf(k, i, l) = (real(l*(l - 1), dp)*rk**(l - 2) - &
629 2._dp*al*real(2*l + 1, dp)*rk**(l) + 4._dp*al*rk**(l + 2))*ear
630 END DO
631 END DO
632 END DO
633
634 CALL set_atom(atom, basis=basis)
635
636 ! optimization defaults
637 atom%optimization%damping = 0.2_dp
638 atom%optimization%eps_scf = 1.e-6_dp
639 atom%optimization%eps_diis = 100._dp
640 atom%optimization%max_iter = 50
641 atom%optimization%n_diis = 5
642
643 nelem = 0
644 ncore = 0
645 ncalc = 0
646 IF (ASSOCIATED(gth_potential)) THEN
647 CALL get_potential(gth_potential, elec_conf=econf)
648 CALL set_pseudo_state(econf, z, ncalc, ncore, nelem)
649 ELSE IF (ASSOCIATED(sgp_potential)) THEN
650 CALL get_potential(sgp_potential, elec_conf=econf)
651 CALL set_pseudo_state(econf, z, ncalc, ncore, nelem)
652 ELSE
653 DO l = 0, min(lmat, ubound(ptable(z)%e_conv, 1))
654 ll = 2*(2*l + 1)
655 nn = ptable(z)%e_conv(l)
656 ii = 0
657 DO
658 ii = ii + 1
659 IF (nn <= ll) THEN
660 nelem(l, ii) = nn
661 EXIT
662 ELSE
663 nelem(l, ii) = ll
664 nn = nn - ll
665 END IF
666 END DO
667 END DO
668 ncalc = nelem - ncore
669 END IF
670
671 IF (qs_kind%ghost .OR. qs_kind%floating) THEN
672 nelem = 0
673 ncore = 0
674 ncalc = 0
675 END IF
676
677 ALLOCATE (atom%state)
678
679 atom%state%core = 0._dp
680 atom%state%core(0:lmat, 1:7) = real(ncore(0:lmat, 1:7), dp)
681 atom%state%occ = 0._dp
682 atom%state%occ(0:lmat, 1:7) = real(ncalc(0:lmat, 1:7), dp)
683 atom%state%occupation = 0._dp
684 atom%state%multiplicity = -1
685 DO l = 0, lmat
686 k = 0
687 DO i = 1, 7
688 IF (ncalc(l, i) > 0) THEN
689 k = k + 1
690 atom%state%occupation(l, k) = real(ncalc(l, i), dp)
691 END IF
692 END DO
693 END DO
694
695 atom%state%maxl_occ = get_maxl_occ(atom%state%occupation)
696 atom%state%maxn_occ = get_maxn_occ(atom%state%occupation)
697 atom%state%maxl_calc = atom%state%maxl_occ
698 atom%state%maxn_calc = atom%state%maxn_occ
699
700 ! calculate integrals
701 ! general integrals
702 CALL atom_int_setup(integrals, basis, potential=atom%potential, &
703 eri_coulomb=(atom%coulomb_integral_type == do_analytic), &
704 eri_exchange=(atom%exchange_integral_type == do_analytic))
705 ! potential
706 CALL atom_ppint_setup(integrals, basis, potential=atom%potential)
707 ! relativistic correction terms
708 NULLIFY (integrals%tzora, integrals%hdkh)
709 CALL atom_relint_setup(integrals, basis, atom%relativistic, zcore=real(atom%zcore, dp))
710 CALL set_atom(atom, integrals=integrals)
711
712 NULLIFY (orbitals)
713 mo = maxval(atom%state%maxn_calc)
714 mb = maxval(atom%basis%nbas)
715 CALL create_atom_orbs(orbitals, mb, mo)
716 CALL set_atom(atom, orbitals=orbitals)
717
718 IF (PRESENT(iunit)) THEN
719 CALL calculate_atom(atom, iunit)
720 CALL atom_fit_density(atom, ngto, 0, iunit, results=results)
721 ELSE
722 CALL calculate_atom(atom, -1)
723 CALL atom_fit_density(atom, ngto, 0, -1, results=results)
724 END IF
725
726 xx = results(1)
727 cc = results(2)
728 DO i = 1, ngto
729 density(i, 1) = xx*cc**i
730 density(i, 2) = results(2 + i)
731 END DO
732
733 ! clean up
734 CALL atom_int_release(integrals)
735 CALL atom_ppint_release(integrals)
736 CALL atom_relint_release(integrals)
737 CALL release_atom_basis(basis)
738 CALL release_atom_potential(potential)
739 CALL release_atom_type(atom)
740
741 DEALLOCATE (potential, basis, integrals)
742
743 END SUBROUTINE calculate_atomic_density
744
745! **************************************************************************************************
746!> \brief ...
747!> \param atomic_kind ...
748!> \param qs_kind ...
749!> \param rel_control ...
750!> \param rtmat ...
751! **************************************************************************************************
752 SUBROUTINE calculate_atomic_relkin(atomic_kind, qs_kind, rel_control, rtmat)
753 TYPE(atomic_kind_type), INTENT(IN) :: atomic_kind
754 TYPE(qs_kind_type), INTENT(IN) :: qs_kind
755 TYPE(rel_control_type), POINTER :: rel_control
756 REAL(kind=dp), DIMENSION(:, :), POINTER :: rtmat
757
758 INTEGER :: i, ii, ipgf, j, k, k1, k2, l, ll, m, n, &
759 ngp, nj, nn, nr, ns, nset, nsgf, &
760 quadtype, relativistic, z
761 INTEGER, DIMENSION(0:lmat, 10) :: ncalc, ncore, nelem
762 INTEGER, DIMENSION(0:lmat, 100) :: set_index, shell_index
763 INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, nshell
764 INTEGER, DIMENSION(:, :), POINTER :: first_sgf, last_sgf, ls
765 REAL(kind=dp) :: al, alpha, ear, prefac, rk, zeff
766 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: omat
767 REAL(kind=dp), DIMENSION(:, :), POINTER :: zet
768 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: gcc
769 TYPE(all_potential_type), POINTER :: all_potential
770 TYPE(atom_basis_type), POINTER :: basis
771 TYPE(atom_integrals), POINTER :: integrals
772 TYPE(atom_potential_type), POINTER :: potential
773 TYPE(atom_type), POINTER :: atom
774 TYPE(grid_atom_type), POINTER :: grid
775 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
776
777 IF (rel_control%rel_method == rel_none) RETURN
778
779 NULLIFY (all_potential, orb_basis_set)
780 CALL get_qs_kind(qs_kind, basis_set=orb_basis_set, all_potential=all_potential)
781
782 cpassert(ASSOCIATED(orb_basis_set))
783
784 IF (ASSOCIATED(all_potential)) THEN
785 ! only all electron atoms will get the relativistic correction
786
787 CALL get_atomic_kind(atomic_kind, z=z)
788 CALL get_qs_kind(qs_kind, zeff=zeff)
789 NULLIFY (atom)
790 CALL create_atom_type(atom)
791 NULLIFY (atom%xc_section)
792 NULLIFY (atom%orbitals)
793 atom%z = z
794 alpha = sqrt(all_potential%alpha_core_charge)
795
796 ! set the method flag
797 SELECT CASE (rel_control%rel_method)
798 CASE DEFAULT
799 cpabort("")
800 CASE (rel_dkh)
801 SELECT CASE (rel_control%rel_DKH_order)
802 CASE DEFAULT
803 cpabort("")
804 CASE (0)
805 relativistic = do_dkh0_atom
806 CASE (1)
807 relativistic = do_dkh1_atom
808 CASE (2)
809 relativistic = do_dkh2_atom
810 CASE (3)
811 relativistic = do_dkh3_atom
812 END SELECT
813 CASE (rel_zora)
814 SELECT CASE (rel_control%rel_zora_type)
815 CASE DEFAULT
816 cpabort("")
817 CASE (rel_zora_full)
818 cpabort("")
819 CASE (rel_zora_mp)
820 relativistic = do_zoramp_atom
821 CASE (rel_sczora_mp)
822 relativistic = do_sczoramp_atom
823 END SELECT
824 END SELECT
825
826 CALL set_atom(atom, &
827 pp_calc=.false., &
828 method_type=do_rks_atom, &
829 relativistic=relativistic, &
830 coulomb_integral_type=do_numeric, &
831 exchange_integral_type=do_numeric)
832
833 ALLOCATE (potential, basis, integrals)
834
835 potential%ppot_type = no_pseudo
836 CALL set_atom(atom, zcore=z, potential=potential)
837
838 CALL get_gto_basis_set(orb_basis_set, &
839 nset=nset, nshell=nshell, npgf=npgf, lmin=lmin, lmax=lmax, l=ls, nsgf=nsgf, zet=zet, gcc=gcc, &
840 first_sgf=first_sgf, last_sgf=last_sgf)
841
842 NULLIFY (grid)
843 ngp = 400
844 quadtype = do_gapw_log
845 CALL allocate_grid_atom(grid)
846 CALL create_grid_atom(grid, ngp, 1, 1, 0, quadtype)
847 grid%nr = ngp
848 basis%grid => grid
849
850 NULLIFY (basis%am, basis%cm, basis%as, basis%ns, basis%bf, basis%dbf, basis%ddbf)
851 basis%basis_type = cgto_basis
852 basis%eps_eig = 1.e-12_dp
853
854 ! fill in the basis data structures
855 set_index = 0
856 shell_index = 0
857 basis%nprim = 0
858 basis%nbas = 0
859 DO i = 1, nset
860 DO j = lmin(i), min(lmax(i), lmat)
861 basis%nprim(j) = basis%nprim(j) + npgf(i)
862 END DO
863 DO j = 1, nshell(i)
864 l = ls(j, i)
865 IF (l <= lmat) THEN
866 basis%nbas(l) = basis%nbas(l) + 1
867 k = basis%nbas(l)
868 cpassert(k <= 100)
869 set_index(l, k) = i
870 shell_index(l, k) = j
871 END IF
872 END DO
873 END DO
874
875 nj = maxval(basis%nprim)
876 ns = maxval(basis%nbas)
877 ALLOCATE (basis%am(nj, 0:lmat))
878 basis%am = 0._dp
879 ALLOCATE (basis%cm(nj, ns, 0:lmat))
880 basis%cm = 0._dp
881 DO j = 0, lmat
882 nj = 0
883 ns = 0
884 DO i = 1, nset
885 IF (j >= lmin(i) .AND. j <= lmax(i)) THEN
886 DO ipgf = 1, npgf(i)
887 basis%am(nj + ipgf, j) = zet(ipgf, i)
888 END DO
889 DO ii = 1, nshell(i)
890 IF (ls(ii, i) == j) THEN
891 ns = ns + 1
892 DO ipgf = 1, npgf(i)
893 basis%cm(nj + ipgf, ns, j) = gcc(ipgf, ii, i)
894 END DO
895 END IF
896 END DO
897 nj = nj + npgf(i)
898 END IF
899 END DO
900 END DO
901
902 ! Normalization as used in the atomic code
903 ! We have to undo the Quickstep normalization
904 DO j = 0, lmat
905 prefac = 2.0_dp*sqrt(pi/dfac(2*j + 1))
906 DO ipgf = 1, basis%nprim(j)
907 DO ii = 1, basis%nbas(j)
908 basis%cm(ipgf, ii, j) = prefac*basis%cm(ipgf, ii, j)
909 END DO
910 END DO
911 END DO
912
913 ! initialize basis function on a radial grid
914 nr = basis%grid%nr
915 m = maxval(basis%nbas)
916 ALLOCATE (basis%bf(nr, m, 0:lmat))
917 ALLOCATE (basis%dbf(nr, m, 0:lmat))
918 ALLOCATE (basis%ddbf(nr, m, 0:lmat))
919
920 basis%bf = 0._dp
921 basis%dbf = 0._dp
922 basis%ddbf = 0._dp
923 DO l = 0, lmat
924 DO i = 1, basis%nprim(l)
925 al = basis%am(i, l)
926 DO k = 1, nr
927 rk = basis%grid%rad(k)
928 ear = exp(-al*basis%grid%rad(k)**2)
929 DO j = 1, basis%nbas(l)
930 basis%bf(k, j, l) = basis%bf(k, j, l) + rk**l*ear*basis%cm(i, j, l)
931 basis%dbf(k, j, l) = basis%dbf(k, j, l) &
932 + (real(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear*basis%cm(i, j, l)
933 basis%ddbf(k, j, l) = basis%ddbf(k, j, l) + &
934 (real(l*(l - 1), dp)*rk**(l - 2) - 2._dp*al*real(2*l + 1, dp)* &
935 rk**(l) + 4._dp*al*rk**(l + 2))*ear*basis%cm(i, j, l)
936 END DO
937 END DO
938 END DO
939 END DO
940
941 CALL set_atom(atom, basis=basis)
942
943 ! optimization defaults
944 atom%optimization%damping = 0.2_dp
945 atom%optimization%eps_scf = 1.e-6_dp
946 atom%optimization%eps_diis = 100._dp
947 atom%optimization%max_iter = 50
948 atom%optimization%n_diis = 5
949
950 ! electronic state
951 nelem = 0
952 ncore = 0
953 ncalc = 0
954 DO l = 0, min(lmat, ubound(ptable(z)%e_conv, 1))
955 ll = 2*(2*l + 1)
956 nn = ptable(z)%e_conv(l)
957 ii = 0
958 DO
959 ii = ii + 1
960 IF (nn <= ll) THEN
961 nelem(l, ii) = nn
962 EXIT
963 ELSE
964 nelem(l, ii) = ll
965 nn = nn - ll
966 END IF
967 END DO
968 END DO
969 ncalc = nelem - ncore
970
971 IF (qs_kind%ghost .OR. qs_kind%floating) THEN
972 nelem = 0
973 ncore = 0
974 ncalc = 0
975 END IF
976
977 ALLOCATE (atom%state)
978
979 atom%state%core = 0._dp
980 atom%state%core(0:lmat, 1:7) = real(ncore(0:lmat, 1:7), dp)
981 atom%state%occ = 0._dp
982 atom%state%occ(0:lmat, 1:7) = real(ncalc(0:lmat, 1:7), dp)
983 atom%state%occupation = 0._dp
984 atom%state%multiplicity = -1
985 DO l = 0, lmat
986 k = 0
987 DO i = 1, 7
988 IF (ncalc(l, i) > 0) THEN
989 k = k + 1
990 atom%state%occupation(l, k) = real(ncalc(l, i), dp)
991 END IF
992 END DO
993 END DO
994
995 atom%state%maxl_occ = get_maxl_occ(atom%state%occupation)
996 atom%state%maxn_occ = get_maxn_occ(atom%state%occupation)
997 atom%state%maxl_calc = atom%state%maxl_occ
998 atom%state%maxn_calc = atom%state%maxn_occ
999
1000 ! calculate integrals
1001 ! general integrals
1002 CALL atom_int_setup(integrals, basis)
1003 ! potential
1004 CALL atom_ppint_setup(integrals, basis, potential=atom%potential)
1005 ! relativistic correction terms
1006 NULLIFY (integrals%tzora, integrals%hdkh)
1007 CALL atom_relint_setup(integrals, basis, atom%relativistic, zcore=real(atom%zcore, dp), &
1008 alpha=alpha)
1009 CALL set_atom(atom, integrals=integrals)
1010
1011 ! for DKH we need erfc integrals to correct non-relativistic
1012 integrals%core = 0.0_dp
1013 DO l = 0, lmat
1014 n = integrals%n(l)
1015 m = basis%nprim(l)
1016 ALLOCATE (omat(m, m))
1017
1018 CALL sg_erfc(omat(1:m, 1:m), l, alpha, basis%am(1:m, l), basis%am(1:m, l))
1019 integrals%core(1:n, 1:n, l) = matmul(transpose(basis%cm(1:m, 1:n, l)), &
1020 matmul(omat(1:m, 1:m), basis%cm(1:m, 1:n, l)))
1021
1022 DEALLOCATE (omat)
1023 END DO
1024
1025 ! recover relativistic kinetic matrix in CP2K/GPW order and normalization
1026 IF (ASSOCIATED(rtmat)) THEN
1027 DEALLOCATE (rtmat)
1028 END IF
1029 ALLOCATE (rtmat(nsgf, nsgf))
1030 rtmat = 0._dp
1031 DO l = 0, lmat
1032 ll = 2*l
1033 DO k1 = 1, basis%nbas(l)
1034 DO k2 = 1, basis%nbas(l)
1035 i = first_sgf(shell_index(l, k1), set_index(l, k1))
1036 j = first_sgf(shell_index(l, k2), set_index(l, k2))
1037 SELECT CASE (atom%relativistic)
1038 CASE DEFAULT
1039 cpabort("")
1040 CASE (do_zoramp_atom, do_sczoramp_atom)
1041 DO m = 0, ll
1042 rtmat(i + m, j + m) = integrals%tzora(k1, k2, l)
1043 END DO
1044 CASE (do_dkh0_atom, do_dkh1_atom, do_dkh2_atom, do_dkh3_atom)
1045 DO m = 0, ll
1046 rtmat(i + m, j + m) = integrals%hdkh(k1, k2, l) - integrals%kin(k1, k2, l) + &
1047 atom%zcore*integrals%core(k1, k2, l)
1048 END DO
1049 END SELECT
1050 END DO
1051 END DO
1052 END DO
1053 DO k1 = 1, nsgf
1054 DO k2 = k1, nsgf
1055 rtmat(k1, k2) = 0.5_dp*(rtmat(k1, k2) + rtmat(k2, k1))
1056 rtmat(k2, k1) = rtmat(k1, k2)
1057 END DO
1058 END DO
1059
1060 ! clean up
1061 CALL atom_int_release(integrals)
1062 CALL atom_ppint_release(integrals)
1063 CALL atom_relint_release(integrals)
1064 CALL release_atom_basis(basis)
1065 CALL release_atom_potential(potential)
1066 CALL release_atom_type(atom)
1067
1068 DEALLOCATE (potential, basis, integrals)
1069
1070 ELSE
1071
1072 IF (ASSOCIATED(rtmat)) THEN
1073 DEALLOCATE (rtmat)
1074 END IF
1075 NULLIFY (rtmat)
1076
1077 END IF
1078
1079 END SUBROUTINE calculate_atomic_relkin
1080
1081! **************************************************************************************************
1082!> \brief ...
1083!> \param gth_potential ...
1084!> \param gth_atompot ...
1085! **************************************************************************************************
1086 SUBROUTINE gth_potential_conversion(gth_potential, gth_atompot)
1087 TYPE(gth_potential_type), POINTER :: gth_potential
1088 TYPE(atom_gthpot_type) :: gth_atompot
1089
1090 INTEGER :: i, j, l, lm, n, ne, nexp_lpot, nexp_lsd, &
1091 nexp_nlcc
1092 INTEGER, DIMENSION(:), POINTER :: nct_lpot, nct_lsd, nct_nlcc, nppnl, &
1093 ppeconf
1094 LOGICAL :: lpot_present, lsd_present, nlcc_present
1095 REAL(kind=dp) :: ac, zeff
1096 REAL(kind=dp), DIMENSION(:), POINTER :: alpha_lpot, alpha_lsd, alpha_nlcc, ap, ce
1097 REAL(kind=dp), DIMENSION(:, :), POINTER :: cval_lpot, cval_lsd, cval_nlcc
1098 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: hp
1099
1100 CALL get_potential(gth_potential, &
1101 zeff=zeff, &
1102 elec_conf=ppeconf, &
1103 alpha_core_charge=ac, &
1104 nexp_ppl=ne, &
1105 cexp_ppl=ce, &
1106 lppnl=lm, &
1107 nprj_ppnl=nppnl, &
1108 alpha_ppnl=ap, &
1109 hprj_ppnl=hp)
1110
1111 gth_atompot%zion = zeff
1112 gth_atompot%rc = sqrt(0.5_dp/ac)
1113 gth_atompot%ncl = ne
1114 gth_atompot%cl(:) = 0._dp
1115 IF (ac > 0._dp) THEN
1116 DO i = 1, ne
1117 gth_atompot%cl(i) = ce(i)/(2._dp*ac)**(i - 1)
1118 END DO
1119 END IF
1120 !extended type
1121 gth_atompot%lpotextended = .false.
1122 gth_atompot%lsdpot = .false.
1123 gth_atompot%nlcc = .false.
1124 gth_atompot%nexp_lpot = 0
1125 gth_atompot%nexp_lsd = 0
1126 gth_atompot%nexp_nlcc = 0
1127 CALL get_potential(gth_potential, &
1128 lpot_present=lpot_present, &
1129 lsd_present=lsd_present, &
1130 nlcc_present=nlcc_present)
1131 IF (lpot_present) THEN
1132 CALL get_potential(gth_potential, &
1133 nexp_lpot=nexp_lpot, &
1134 alpha_lpot=alpha_lpot, &
1135 nct_lpot=nct_lpot, &
1136 cval_lpot=cval_lpot)
1137 gth_atompot%lpotextended = .true.
1138 gth_atompot%nexp_lpot = nexp_lpot
1139 gth_atompot%alpha_lpot(1:nexp_lpot) = sqrt(0.5_dp/alpha_lpot(1:nexp_lpot))
1140 gth_atompot%nct_lpot(1:nexp_lpot) = nct_lpot(1:nexp_lpot)
1141 DO j = 1, nexp_lpot
1142 ac = alpha_lpot(j)
1143 DO i = 1, 4
1144 gth_atompot%cval_lpot(i, j) = cval_lpot(i, j)/(2._dp*ac)**(i - 1)
1145 END DO
1146 END DO
1147 END IF
1148 IF (lsd_present) THEN
1149 CALL get_potential(gth_potential, &
1150 nexp_lsd=nexp_lsd, &
1151 alpha_lsd=alpha_lsd, &
1152 nct_lsd=nct_lsd, &
1153 cval_lsd=cval_lsd)
1154 gth_atompot%lsdpot = .true.
1155 gth_atompot%nexp_lsd = nexp_lsd
1156 gth_atompot%alpha_lsd(1:nexp_lsd) = sqrt(0.5_dp/alpha_lsd(1:nexp_lsd))
1157 gth_atompot%nct_lsd(1:nexp_lsd) = nct_lsd(1:nexp_lsd)
1158 DO j = 1, nexp_lpot
1159 ac = alpha_lsd(j)
1160 DO i = 1, 4
1161 gth_atompot%cval_lsd(i, j) = cval_lsd(i, j)/(2._dp*ac)**(i - 1)
1162 END DO
1163 END DO
1164 END IF
1165
1166 ! nonlocal part
1167 gth_atompot%nl(:) = 0
1168 gth_atompot%rcnl(:) = 0._dp
1169 gth_atompot%hnl(:, :, :) = 0._dp
1170 DO l = 0, lm
1171 n = nppnl(l)
1172 gth_atompot%nl(l) = n
1173 gth_atompot%rcnl(l) = sqrt(0.5_dp/ap(l))
1174 gth_atompot%hnl(1:n, 1:n, l) = hp(1:n, 1:n, l)
1175 END DO
1176
1177 IF (nlcc_present) THEN
1178 CALL get_potential(gth_potential, &
1179 nexp_nlcc=nexp_nlcc, &
1180 alpha_nlcc=alpha_nlcc, &
1181 nct_nlcc=nct_nlcc, &
1182 cval_nlcc=cval_nlcc)
1183 gth_atompot%nlcc = .true.
1184 gth_atompot%nexp_nlcc = nexp_nlcc
1185 gth_atompot%alpha_nlcc(1:nexp_nlcc) = alpha_nlcc(1:nexp_nlcc)
1186 gth_atompot%nct_nlcc(1:nexp_nlcc) = nct_nlcc(1:nexp_nlcc)
1187 gth_atompot%cval_nlcc(1:4, 1:nexp_nlcc) = cval_nlcc(1:4, 1:nexp_nlcc)
1188 END IF
1189
1190 END SUBROUTINE gth_potential_conversion
1191
1192! **************************************************************************************************
1193!> \brief ...
1194!> \param sgp_potential ...
1195!> \param sgp_atompot ...
1196! **************************************************************************************************
1197 SUBROUTINE sgp_potential_conversion(sgp_potential, sgp_atompot)
1198 TYPE(sgp_potential_type), POINTER :: sgp_potential
1199 TYPE(atom_sgppot_type) :: sgp_atompot
1200
1201 INTEGER :: lm, n
1202 INTEGER, DIMENSION(:), POINTER :: ppeconf
1203 LOGICAL :: nlcc_present
1204 REAL(kind=dp) :: ac, zeff
1205 REAL(kind=dp), DIMENSION(:), POINTER :: ap, ce
1206 REAL(kind=dp), DIMENSION(:, :), POINTER :: hhp
1207 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: ccp
1208
1209 CALL get_potential(sgp_potential, &
1210 name=sgp_atompot%pname, &
1211 zeff=zeff, &
1212 elec_conf=ppeconf, &
1213 alpha_core_charge=ac)
1214 sgp_atompot%zion = zeff
1215 sgp_atompot%ac_local = ac
1216 sgp_atompot%econf(0:3) = ppeconf(0:3)
1217 CALL get_potential(sgp_potential, lmax=lm, &
1218 is_nonlocal=sgp_atompot%is_nonlocal, &
1219 n_nonlocal=n, a_nonlocal=ap, h_nonlocal=hhp, c_nonlocal=ccp)
1220 ! nonlocal
1221 sgp_atompot%has_nonlocal = any(sgp_atompot%is_nonlocal)
1222 sgp_atompot%lmax = lm
1223 IF (sgp_atompot%has_nonlocal) THEN
1224 cpassert(n <= SIZE(sgp_atompot%a_nonlocal))
1225 sgp_atompot%n_nonlocal = n
1226 sgp_atompot%a_nonlocal(1:n) = ap(1:n)
1227 sgp_atompot%h_nonlocal(1:n, 0:lm) = hhp(1:n, 0:lm)
1228 sgp_atompot%c_nonlocal(1:n, 1:n, 0:lm) = ccp(1:n, 1:n, 0:lm)
1229 END IF
1230 ! local
1231 CALL get_potential(sgp_potential, n_local=n, a_local=ap, c_local=ce)
1232 cpassert(n <= SIZE(sgp_atompot%a_local))
1233 sgp_atompot%n_local = n
1234 sgp_atompot%a_local(1:n) = ap(1:n)
1235 sgp_atompot%c_local(1:n) = ce(1:n)
1236 ! NLCC
1237 CALL get_potential(sgp_potential, has_nlcc=nlcc_present, &
1238 n_nlcc=n, a_nlcc=ap, c_nlcc=ce)
1239 IF (nlcc_present) THEN
1240 sgp_atompot%has_nlcc = .true.
1241 cpassert(n <= SIZE(sgp_atompot%a_nlcc))
1242 sgp_atompot%n_nlcc = n
1243 sgp_atompot%a_nlcc(1:n) = ap(1:n)
1244 sgp_atompot%c_nlcc(1:n) = ce(1:n)
1245 ELSE
1246 sgp_atompot%has_nlcc = .false.
1247 END IF
1248
1249 END SUBROUTINE sgp_potential_conversion
1250
1251! **************************************************************************************************
1252!> \brief ...
1253!> \param sgp_potential ...
1254!> \param ecp_atompot ...
1255! **************************************************************************************************
1256 SUBROUTINE ecp_potential_conversion(sgp_potential, ecp_atompot)
1257 TYPE(sgp_potential_type), POINTER :: sgp_potential
1258 TYPE(atom_ecppot_type) :: ecp_atompot
1259
1260 INTEGER, DIMENSION(:), POINTER :: ppeconf
1261 LOGICAL :: ecp_local, ecp_semi_local
1262 REAL(kind=dp) :: zeff
1263
1264 CALL get_potential(sgp_potential, ecp_local=ecp_local, ecp_semi_local=ecp_semi_local)
1265 cpassert(ecp_semi_local .AND. ecp_local)
1266 CALL get_potential(sgp_potential, &
1267 name=ecp_atompot%pname, &
1268 zeff=zeff, &
1269 elec_conf=ppeconf)
1270 ecp_atompot%zion = zeff
1271 ecp_atompot%econf(0:3) = ppeconf(0:3)
1272 CALL get_potential(sgp_potential, lmax=ecp_atompot%lmax)
1273 ! local
1274 CALL get_potential(sgp_potential, nloc=ecp_atompot%nloc, nrloc=ecp_atompot%nrloc, &
1275 aloc=ecp_atompot%aloc, bloc=ecp_atompot%bloc)
1276 ! nonlocal
1277 CALL get_potential(sgp_potential, npot=ecp_atompot%npot, nrpot=ecp_atompot%nrpot, &
1278 apot=ecp_atompot%apot, bpot=ecp_atompot%bpot)
1279
1280 END SUBROUTINE ecp_potential_conversion
1281! **************************************************************************************************
1282
1283END MODULE atom_kind_orbitals
external_potential_types::get_potential
Definition external_potential_types.F:276
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, powell_section, results)
Fit the atomic electron density using a geometrical Gaussian basis set.
Definition atom_fit.F:77
atom_kind_orbitals
calculate the orbitals for a given atomic kind type
Definition atom_kind_orbitals.F:11
atom_kind_orbitals::calculate_atomic_density
subroutine, public calculate_atomic_density(density, atomic_kind, qs_kind, ngto, iunit, allelectron, confine)
...
Definition atom_kind_orbitals.F:479
atom_kind_orbitals::calculate_atomic_relkin
subroutine, public calculate_atomic_relkin(atomic_kind, qs_kind, rel_control, rtmat)
...
Definition atom_kind_orbitals.F:753
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::gth_potential_conversion
subroutine, public gth_potential_conversion(gth_potential, gth_atompot)
...
Definition atom_kind_orbitals.F:1087
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:1145
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:64
atom_operators::atom_relint_setup
subroutine, public atom_relint_setup(integrals, basis, reltyp, zcore, alpha)
...
Definition atom_operators.F:802
atom_operators::atom_relint_release
subroutine, public atom_relint_release(integrals)
Release memory allocated for atomic integrals (relativistic effects).
Definition atom_operators.F:1166
atom_operators::atom_ppint_setup
subroutine, public atom_ppint_setup(integrals, basis, potential)
...
Definition atom_operators.F:337
atom_operators::atom_int_release
subroutine, public atom_int_release(integrals)
Release memory allocated for atomic integrals (valence electrons).
Definition atom_operators.F:1108
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:944
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:968
atom_types::release_atom_potential
subroutine, public release_atom_potential(potential)
...
Definition atom_types.F:2521
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:1009
atom_types::release_atom_basis
subroutine, public release_atom_basis(basis)
...
Definition atom_types.F:910
atom_types::create_atom_orbs
subroutine, public create_atom_orbs(orbs, mbas, mo)
...
Definition atom_types.F:1055
atom_types::clementi_geobas
subroutine, public clementi_geobas(zval, cval, aval, ngto, ival)
...
Definition atom_types.F:1428
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:302
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:367
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:282
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)
...
Definition basis_set_types.F:615
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:606
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:606
input_constants::rel_zora
integer, parameter, public rel_zora
Definition input_constants.F:603
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:603
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:225
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:603
input_constants::rel_sczora_mp
integer, parameter, public rel_sczora_mp
Definition input_constants.F:606
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:3285
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:3499
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, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, 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_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:451
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:3152
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:130
atom_types::atom_gthpot_type
Provides all information about a pseudopotential.
Definition atom_types.F:98
atom_types::atom_integrals
Definition atom_types.F:213
atom_types::atom_orbitals
Holds atomic orbitals and energies.
Definition atom_types.F:234
atom_types::atom_potential_type
Definition atom_types.F:168
atom_types::atom_sgppot_type
Definition atom_types.F:146
atom_types::atom_type
Provides all information about an atomic kind.
Definition atom_types.F:290
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:71
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:172
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126
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:171
rel_control_types::rel_control_type
contains the parameters needed by a relativistic calculation
Definition rel_control_types.F:55