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