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