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