atom_kind_orbitals.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief calculate the orbitals for a given atomic kind type
! **************************************************************************************************
MODULE atom_kind_orbitals
   USE ai_onecenter,                    ONLY: sg_erfc
   USE atom_electronic_structure,       ONLY: calculate_atom
   USE atom_fit,                        ONLY: atom_fit_density
   USE atom_operators,                  ONLY: atom_int_release,&
                                              atom_int_setup,&
                                              atom_ppint_release,&
                                              atom_ppint_setup,&
                                              atom_relint_release,&
                                              atom_relint_setup
   USE atom_set_basis,                  ONLY: set_kind_basis_atomic
   USE atom_types,                      ONLY: &
        cgto_basis, clementi_geobas, gto_basis, atom_basis_type, atom_ecppot_type, &
        atom_gthpot_type, atom_integrals, atom_orbitals, atom_potential_type, atom_sgppot_type, &
        atom_type, create_atom_orbs, create_atom_type, lmat, release_atom_basis, &
        release_atom_potential, release_atom_type, set_atom
   USE atom_utils,                      ONLY: atom_density,&
                                              get_maxl_occ,&
                                              get_maxn_occ
   USE atomic_kind_types,               ONLY: atomic_kind_type,&
                                              get_atomic_kind
   USE basis_set_types,                 ONLY: get_gto_basis_set,&
                                              gto_basis_set_type
   USE external_potential_types,        ONLY: all_potential_type,&
                                              get_potential,&
                                              gth_potential_type,&
                                              sgp_potential_type
   USE input_constants,                 ONLY: &
        barrier_conf, do_analytic, do_dkh0_atom, do_dkh1_atom, do_dkh2_atom, do_dkh3_atom, &
        do_gapw_log, do_nonrel_atom, do_numeric, do_rks_atom, do_sczoramp_atom, do_uks_atom, &
        do_zoramp_atom, ecp_pseudo, gth_pseudo, no_pseudo, poly_conf, rel_dkh, rel_none, &
        rel_sczora_mp, rel_zora, rel_zora_full, rel_zora_mp, sgp_pseudo, upf_pseudo
   USE input_section_types,             ONLY: section_vals_type
   USE kinds,                           ONLY: dp
   USE mathconstants,                   ONLY: dfac,&
                                              fourpi,&
                                              pi
   USE periodic_table,                  ONLY: ptable
   USE physcon,                         ONLY: bohr
   USE qs_grid_atom,                    ONLY: allocate_grid_atom,&
                                              create_grid_atom,&
                                              grid_atom_type
   USE qs_kind_types,                   ONLY: get_qs_kind,&
                                              init_atom_electronic_state,&
                                              qs_kind_type,&
                                              set_pseudo_state
   USE rel_control_types,               ONLY: rel_control_type
#include "./base/base_uses.f90"
 
   IMPLICIT NONE
 
   PRIVATE
 
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_kind_orbitals'
 
   PUBLIC :: calculate_atomic_orbitals, calculate_atomic_density, &
             calculate_atomic_relkin, gth_potential_conversion, sgp_potential_conversion
 
! **************************************************************************************************
 
CONTAINS
 
! **************************************************************************************************
!> \brief ...
!> \param atomic_kind ...
!> \param qs_kind ...
!> \param agrid ...
!> \param iunit ...
!> \param pmat ...
!> \param fmat ...
!> \param density ...
!> \param wavefunction ...
!> \param wfninfo ...
!> \param confine ...
!> \param xc_section ...
!> \param nocc ...
!> \param which_l ...
!> \param which_n ...
!> \param proj_shell_charge ...
!> \param ao_coef ...
! **************************************************************************************************
   SUBROUTINE 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)
      TYPE(atomic_kind_type), INTENT(IN)                 :: atomic_kind
      TYPE(qs_kind_type), INTENT(IN)                     :: qs_kind
      TYPE(grid_atom_type), OPTIONAL                     :: agrid
      INTEGER, INTENT(IN), OPTIONAL                      :: iunit
      REAL(kind=dp), DIMENSION(:, :, :), OPTIONAL, &
         POINTER                                         :: pmat, fmat
      REAL(kind=dp), DIMENSION(:), OPTIONAL, POINTER     :: density
      REAL(kind=dp), DIMENSION(:, :), OPTIONAL, POINTER  :: wavefunction, wfninfo
      LOGICAL, INTENT(IN), OPTIONAL                      :: confine
      TYPE(section_vals_type), OPTIONAL, POINTER         :: xc_section
      INTEGER, DIMENSION(:), OPTIONAL                    :: nocc
      INTEGER, INTENT(IN), OPTIONAL                      :: which_l, which_n
      REAL(kind=dp), INTENT(IN), OPTIONAL                :: proj_shell_charge(4)
      REAL(kind=dp), DIMENSION(:), INTENT(INOUT), &
         OPTIONAL, POINTER                               :: ao_coef
 
      INTEGER                                            :: i, ii, j, k, k1, k2, l, ll, m, mb, mo, &
                                                            nr, nset, nsgf, projector_index, &
                                                            valence_n, z
      INTEGER, DIMENSION(0:lmat)                         :: nbb
      INTEGER, DIMENSION(0:lmat, 10)                     :: ncalc, ncore, nelem
      INTEGER, DIMENSION(0:lmat, 100)                    :: set_index, shell_index
      INTEGER, DIMENSION(:), POINTER                     :: nshell
      INTEGER, DIMENSION(:, :), POINTER                  :: first_sgf, ls
      LOGICAL                                            :: ecp_semi_local, ghost, has_pp, uks
      REAL(kind=dp)                                      :: ok, qleft, scal, take, zeff
      REAL(kind=dp), DIMENSION(0:lmat, 10)               :: occupation_neutral
      REAL(kind=dp), DIMENSION(0:lmat, 10, 2)            :: edelta
      TYPE(all_potential_type), POINTER                  :: all_potential
      TYPE(atom_basis_type), POINTER                     :: basis
      TYPE(atom_integrals), POINTER                      :: integrals
      TYPE(atom_orbitals), POINTER                       :: orbitals
      TYPE(atom_potential_type), POINTER                 :: potential
      TYPE(atom_type), POINTER                           :: atom
      TYPE(gth_potential_type), POINTER                  :: gth_potential
      TYPE(gto_basis_set_type), POINTER                  :: orb_basis_set
      TYPE(sgp_potential_type), POINTER                  :: sgp_potential
 
      IF (PRESENT(ao_coef) .OR. PRESENT(which_l) .OR. &
          PRESENT(which_n) .OR. PRESENT(proj_shell_charge)) THEN
         IF (.NOT. PRESENT(ao_coef) .OR. .NOT. PRESENT(which_l) .OR. &
             .NOT. PRESENT(which_n) .OR. .NOT. PRESENT(proj_shell_charge)) THEN
            cpabort("ao_coef, which_l, which_n, and proj_shell_charge must be provided together")
         END IF
      END IF
 
      NULLIFY (atom)
      CALL create_atom_type(atom)
 
      IF (PRESENT(xc_section)) THEN
         atom%xc_section => xc_section
      ELSE
         NULLIFY (atom%xc_section)
      END IF
 
      CALL get_atomic_kind(atomic_kind, z=z)
      NULLIFY (all_potential, gth_potential, sgp_potential, orb_basis_set)
      CALL get_qs_kind(qs_kind, zeff=zeff, &
                       basis_set=orb_basis_set, &
                       ghost=ghost, &
                       all_potential=all_potential, &
                       gth_potential=gth_potential, &
                       sgp_potential=sgp_potential)
 
      has_pp = ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)
 
      atom%z = z
      CALL set_atom(atom, &
                    pp_calc=has_pp, &
                    do_zmp=.false., &
                    doread=.false., &
                    read_vxc=.false., &
                    relativistic=do_nonrel_atom, &
                    coulomb_integral_type=do_numeric, &
                    exchange_integral_type=do_numeric)
 
      ALLOCATE (potential, integrals)
 
      IF (PRESENT(confine)) THEN
         potential%confinement = confine
      ELSE
         IF (ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)) THEN
            potential%confinement = .true.
         ELSE
            potential%confinement = .false.
         END IF
      END IF
      potential%conf_type = poly_conf
      potential%acon = 0.1_dp
      potential%rcon = 2.0_dp*ptable(z)%vdw_radius*bohr
      potential%scon = 2.0_dp
 
      IF (ASSOCIATED(gth_potential)) THEN
         potential%ppot_type = gth_pseudo
         CALL get_potential(gth_potential, zeff=zeff)
         CALL gth_potential_conversion(gth_potential, potential%gth_pot)
         CALL set_atom(atom, zcore=nint(zeff), potential=potential)
      ELSE IF (ASSOCIATED(sgp_potential)) THEN
         CALL get_potential(sgp_potential, ecp_semi_local=ecp_semi_local)
         IF (ALLOCATED(sgp_potential%upf_data)) THEN
            potential%ppot_type = upf_pseudo
            potential%upf_pot = sgp_potential%upf_data
         ELSE IF (ecp_semi_local) THEN
            potential%ppot_type = ecp_pseudo
            CALL ecp_potential_conversion(sgp_potential, potential%ecp_pot)
            potential%ecp_pot%symbol = ptable(z)%symbol
         ELSE
            potential%ppot_type = sgp_pseudo
            CALL sgp_potential_conversion(sgp_potential, potential%sgp_pot)
            potential%sgp_pot%symbol = ptable(z)%symbol
         END IF
         CALL get_potential(sgp_potential, zeff=zeff)
         CALL set_atom(atom, zcore=nint(zeff), potential=potential)
      ELSE
         potential%ppot_type = no_pseudo
         CALL set_atom(atom, zcore=z, potential=potential)
      END IF
 
      NULLIFY (basis)
      ALLOCATE (basis)
 
      CALL set_kind_basis_atomic(basis, orb_basis_set, has_pp, agrid)
 
      CALL set_atom(atom, basis=basis)
 
      ! optimization defaults
      atom%optimization%damping = 0.2_dp
      atom%optimization%eps_scf = 1.e-6_dp
      atom%optimization%eps_diis = 100._dp
      atom%optimization%max_iter = 50
      atom%optimization%n_diis = 5
 
      ! set up the electronic state
      CALL init_atom_electronic_state(atomic_kind=atomic_kind, &
                                      qs_kind=qs_kind, &
                                      ncalc=ncalc, &
                                      ncore=ncore, &
                                      nelem=nelem, &
                                      edelta=edelta, &
                                      rks=PRESENT(ao_coef))
 
      ! restricted or unrestricted?
      IF (sum(abs(edelta)) > 0.0_dp) THEN
         uks = .true.
         CALL set_atom(atom, method_type=do_uks_atom)
      ELSE
         uks = .false.
         CALL set_atom(atom, method_type=do_rks_atom)
      END IF
 
      ALLOCATE (atom%state)
 
      atom%state%core = 0._dp
      atom%state%core(0:lmat, 1:7) = real(ncore(0:lmat, 1:7), dp)
      atom%state%occ = 0._dp
      IF (uks) THEN
         atom%state%occ(0:lmat, 1:7) = real(ncalc(0:lmat, 1:7), dp) + &
                                       edelta(0:lmat, 1:7, 1) + edelta(0:lmat, 1:7, 2)
      ELSE
         atom%state%occ(0:lmat, 1:7) = real(ncalc(0:lmat, 1:7), dp)
      END IF
      atom%state%occupation = 0._dp
      DO l = 0, lmat
         k = 0
         DO i = 1, 7
            IF (ncalc(l, i) > 0) THEN
               k = k + 1
               IF (uks) THEN
                  atom%state%occupation(l, k) = real(ncalc(l, i), dp) + &
                                                edelta(l, i, 1) + edelta(l, i, 2)
                  atom%state%occa(l, k) = 0.5_dp*real(ncalc(l, i), dp) + edelta(l, i, 1)
                  atom%state%occb(l, k) = 0.5_dp*real(ncalc(l, i), dp) + edelta(l, i, 2)
               ELSE
                  atom%state%occupation(l, k) = real(ncalc(l, i), dp)
               END IF
            END IF
         END DO
         ok = real(2*l + 1, kind=dp)
         IF (uks) THEN
            DO i = 1, 7
               atom%state%occ(l, i) = min(atom%state%occ(l, i), 2.0_dp*ok)
               atom%state%occa(l, i) = min(atom%state%occa(l, i), ok)
               atom%state%occb(l, i) = min(atom%state%occb(l, i), ok)
               atom%state%occupation(l, i) = atom%state%occa(l, i) + atom%state%occb(l, i)
            END DO
         ELSE
            DO i = 1, 7
               atom%state%occ(l, i) = min(atom%state%occ(l, i), 2.0_dp*ok)
               atom%state%occupation(l, i) = min(atom%state%occupation(l, i), 2.0_dp*ok)
            END DO
         END IF
      END DO
 
      IF (PRESENT(proj_shell_charge)) THEN
         ! Save the neutral/original calculated shell list.
         occupation_neutral(0:lmat, 1:10) = atom%state%occupation(0:lmat, 1:10)
         IF (uks) THEN
            CALL cp_abort(__location__, "proj_shell_charge in"// &
                          " calculate_atomic_orbitals is implemented only for RKS")
         END IF
         DO l = 0, min(3, lmat)
            qleft = proj_shell_charge(l + 1)
            IF (qleft > 0.0_dp) THEN
               ! Remove electrons from the highest shell in this angular momentum channel.
               DO k = 7, 1, -1
                  IF (atom%state%occupation(l, k) > 0.0_dp) THEN
                     take = min(qleft, atom%state%occupation(l, k))
                     atom%state%occupation(l, k) = atom%state%occupation(l, k) - take
                     qleft = qleft - take
                     IF (qleft <= 1.0e-12_dp) EXIT
                  END IF
               END DO
               IF (qleft > 1.0e-10_dp) THEN
                  CALL cp_abort(__location__, "proj_shell_charge removes"// &
                                " more electrons than available in this angular channel")
               END IF
               ! Keep atom%state%occ consistent with atom%state%occupation.
               qleft = proj_shell_charge(l + 1)
               DO k = 7, 1, -1
                  IF (atom%state%occ(l, k) > 0.0_dp) THEN
                     take = min(qleft, atom%state%occ(l, k))
                     atom%state%occ(l, k) = atom%state%occ(l, k) - take
                     qleft = qleft - take
                     IF (qleft <= 1.0e-12_dp) EXIT
                  END IF
               END DO
            ELSE IF (qleft < 0.0_dp) THEN
               CALL cp_abort(__location__, "Negative proj_shell_charge"// &
                             " is not supported in this implementation")
            END IF
         END DO
      END IF
 
      IF (uks) THEN
         atom%state%multiplicity = nint(abs(sum(atom%state%occa - atom%state%occb)) + 1)
      ELSE
         atom%state%multiplicity = -1
      END IF
 
      atom%state%maxl_occ = get_maxl_occ(atom%state%occupation)
      atom%state%maxn_occ = get_maxn_occ(atom%state%occupation)
      IF (PRESENT(proj_shell_charge)) THEN
         atom%state%maxl_calc = get_maxl_occ(occupation_neutral)
         atom%state%maxn_calc = get_maxn_occ(occupation_neutral)
      ELSE
         atom%state%maxl_calc = atom%state%maxl_occ
         atom%state%maxn_calc = atom%state%maxn_occ
      END IF
 
      ! total number of occupied orbitals
      IF (PRESENT(nocc) .AND. ghost) THEN
         nocc = 0
      ELSE IF (PRESENT(nocc)) THEN
         nocc = 0
         DO l = 0, lmat
            DO k = 1, 7
               IF (uks) THEN
                  IF (atom%state%occa(l, k) > 0.0_dp) THEN
                     nocc(1) = nocc(1) + 2*l + 1
                  END IF
                  IF (atom%state%occb(l, k) > 0.0_dp) THEN
                     nocc(2) = nocc(2) + 2*l + 1
                  END IF
               ELSE
                  IF (atom%state%occupation(l, k) > 0.0_dp) THEN
                     nocc(1) = nocc(1) + 2*l + 1
                     nocc(2) = nocc(2) + 2*l + 1
                  END IF
               END IF
            END DO
         END DO
      END IF
 
      ! calculate integrals
      ! general integrals
      CALL atom_int_setup(integrals, basis, potential=atom%potential, &
                          eri_coulomb=(atom%coulomb_integral_type == do_analytic), &
                          eri_exchange=(atom%exchange_integral_type == do_analytic))
      ! potential
      CALL atom_ppint_setup(integrals, basis, potential=atom%potential)
      ! relativistic correction terms
      NULLIFY (integrals%tzora, integrals%hdkh)
      CALL atom_relint_setup(integrals, basis, atom%relativistic, zcore=real(atom%zcore, dp))
      CALL set_atom(atom, integrals=integrals)
 
      NULLIFY (orbitals)
      mo = maxval(atom%state%maxn_calc)
      mb = maxval(atom%basis%nbas)
      CALL create_atom_orbs(orbitals, mb, mo)
      CALL set_atom(atom, orbitals=orbitals)
 
      IF (.NOT. ghost) THEN
         IF (PRESENT(iunit)) THEN
            CALL calculate_atom(atom, iunit)
         ELSE
            CALL calculate_atom(atom, -1)
         END IF
      END IF
      IF (PRESENT(pmat)) THEN
         ! recover density matrix in CP2K/GPW order and normalization
         CALL get_gto_basis_set(orb_basis_set, &
                                nset=nset, nshell=nshell, l=ls, nsgf=nsgf, first_sgf=first_sgf)
         set_index = 0
         shell_index = 0
         nbb = 0
         DO i = 1, nset
            DO j = 1, nshell(i)
               l = ls(j, i)
               IF (l <= lmat) THEN
                  nbb(l) = nbb(l) + 1
                  k = nbb(l)
                  cpassert(k <= 100)
                  set_index(l, k) = i
                  shell_index(l, k) = j
               END IF
            END DO
         END DO
 
         IF (ASSOCIATED(pmat)) THEN
            DEALLOCATE (pmat)
         END IF
         ALLOCATE (pmat(nsgf, nsgf, 2))
         pmat = 0._dp
         IF (.NOT. ghost) THEN
            DO l = 0, lmat
               ll = 2*l
               DO k1 = 1, atom%basis%nbas(l)
                  DO k2 = 1, atom%basis%nbas(l)
                     scal = sqrt(atom%integrals%ovlp(k1, k1, l)*atom%integrals%ovlp(k2, k2, l))/real(2*l + 1, kind=dp)
                     i = first_sgf(shell_index(l, k1), set_index(l, k1))
                     j = first_sgf(shell_index(l, k2), set_index(l, k2))
                     IF (uks) THEN
                        DO m = 0, ll
                           pmat(i + m, j + m, 1) = atom%orbitals%pmata(k1, k2, l)*scal
                           pmat(i + m, j + m, 2) = atom%orbitals%pmatb(k1, k2, l)*scal
                        END DO
                     ELSE
                        DO m = 0, ll
                           pmat(i + m, j + m, 1) = atom%orbitals%pmat(k1, k2, l)*scal
                        END DO
                     END IF
                  END DO
               END DO
            END DO
            IF (uks) THEN
               pmat(:, :, 1) = pmat(:, :, 1) + pmat(:, :, 2)
               pmat(:, :, 2) = pmat(:, :, 1) - 2.0_dp*pmat(:, :, 2)
            END IF
         END IF
      END IF
 
      IF (PRESENT(fmat)) THEN
         ! recover fock matrix in CP2K/GPW order.
         ! Caution: Normalization is not take care of, so it's probably weird.
         CALL get_gto_basis_set(orb_basis_set, &
                                nset=nset, nshell=nshell, l=ls, nsgf=nsgf, first_sgf=first_sgf)
         set_index = 0
         shell_index = 0
         nbb = 0
         DO i = 1, nset
            DO j = 1, nshell(i)
               l = ls(j, i)
               IF (l <= lmat) THEN
                  nbb(l) = nbb(l) + 1
                  k = nbb(l)
                  cpassert(k <= 100)
                  set_index(l, k) = i
                  shell_index(l, k) = j
               END IF
            END DO
         END DO
         IF (uks) cpabort("calculate_atomic_orbitals: only RKS is implemented")
         IF (ASSOCIATED(fmat)) cpabort("fmat already associated")
         IF (.NOT. ASSOCIATED(atom%fmat)) cpabort("atom%fmat not associated")
         ALLOCATE (fmat(nsgf, nsgf, 1))
         fmat = 0.0_dp
         IF (.NOT. ghost) THEN
            DO l = 0, lmat
               ll = 2*l
               DO k1 = 1, atom%basis%nbas(l)
               DO k2 = 1, atom%basis%nbas(l)
                  scal = sqrt(atom%integrals%ovlp(k1, k1, l)*atom%integrals%ovlp(k2, k2, l))
                  i = first_sgf(shell_index(l, k1), set_index(l, k1))
                  j = first_sgf(shell_index(l, k2), set_index(l, k2))
                  DO m = 0, ll
                     fmat(i + m, j + m, 1) = atom%fmat%op(k1, k2, l)/scal
                  END DO
               END DO
               END DO
            END DO
         END IF
      END IF
 
      nr = basis%grid%nr
 
      IF (PRESENT(density)) THEN
         IF (ASSOCIATED(density)) DEALLOCATE (density)
         ALLOCATE (density(nr))
         IF (ghost) THEN
            density = 0.0_dp
         ELSE
            CALL atom_density(density, atom%orbitals%pmat, atom%basis, atom%state%maxl_occ)
         END IF
      END IF
 
      IF (PRESENT(wavefunction)) THEN
         cpassert(PRESENT(wfninfo))
         IF (ASSOCIATED(wavefunction)) DEALLOCATE (wavefunction)
         IF (ASSOCIATED(wfninfo)) DEALLOCATE (wfninfo)
         mo = sum(atom%state%maxn_occ)
         ALLOCATE (wavefunction(nr, mo), wfninfo(2, mo))
         wavefunction = 0.0_dp
         IF (.NOT. ghost) THEN
            ii = 0
            DO l = 0, lmat
               DO i = 1, atom%state%maxn_occ(l)
                  IF (atom%state%occupation(l, i) > 0.0_dp) THEN
                     ii = ii + 1
                     wfninfo(1, ii) = atom%state%occupation(l, i)
                     wfninfo(2, ii) = real(l, dp)
                     DO j = 1, atom%basis%nbas(l)
                        wavefunction(:, ii) = wavefunction(:, ii) + &
                                              atom%orbitals%wfn(j, i, l)*basis%bf(:, j, l)
                     END DO
                  END IF
               END DO
            END DO
            cpassert(mo == ii)
         END IF
      END IF
 
      IF (PRESENT(ao_coef)) THEN
         cpassert(PRESENT(which_l))
         cpassert(PRESENT(which_n))
         cpassert(which_l >= 0 .AND. which_l <= lmat)
         valence_n = count(ncore(which_l, :) > 0) + which_l + 1
         projector_index = which_n - valence_n + 1
         IF (projector_index < 1 .OR. projector_index > SIZE(atom%orbitals%wfn, 2)) THEN
            CALL cp_abort(__location__, &
                          "The requested tensorial atomic projector N is outside "// &
                          "the available orbital range for the selected L.")
         END IF
         IF (ASSOCIATED(ao_coef)) DEALLOCATE (ao_coef)
         ALLOCATE (ao_coef(SIZE(atom%orbitals%wfn(:, 1, which_l))))
         ao_coef(:) = atom%orbitals%wfn(:, projector_index, which_l)
         IF (maxval(abs(ao_coef)) <= 100.0_dp*epsilon(1.0_dp)) THEN
            CALL cp_abort(__location__, &
                          "The requested atomic projector is zero. Check the "// &
                          "specified N and L quantum numbers.")
         END IF
         DO j = 1, SIZE(ao_coef(:))
            ao_coef(j) = ao_coef(j)*sqrt(atom%integrals%ovlp(j, j, which_l))
         END DO
      END IF
 
      ! clean up
      CALL atom_int_release(integrals)
      CALL atom_ppint_release(integrals)
      CALL atom_relint_release(integrals)
      CALL release_atom_basis(basis)
      CALL release_atom_potential(potential)
      CALL release_atom_type(atom)
 
      DEALLOCATE (potential, basis, integrals)
 
   END SUBROUTINE calculate_atomic_orbitals
 
! **************************************************************************************************
!> \brief ...
!> \param density ...
!> \param atomic_kind ...
!> \param qs_kind ...
!> \param ngto ...
!> \param iunit ...
!> \param optbasis ... Default=T, if basis should be optimized, if not basis is given in input (density)
!> \param allelectron ...
!> \param confine ...
! **************************************************************************************************
   SUBROUTINE calculate_atomic_density(density, atomic_kind, qs_kind, ngto, iunit, &
                                       optbasis, allelectron, confine)
      REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT)      :: density
      TYPE(atomic_kind_type), POINTER                    :: atomic_kind
      TYPE(qs_kind_type), POINTER                        :: qs_kind
      INTEGER, INTENT(IN)                                :: ngto
      INTEGER, INTENT(IN), OPTIONAL                      :: iunit
      LOGICAL, INTENT(IN), OPTIONAL                      :: optbasis, allelectron, confine
 
      INTEGER, PARAMETER                                 :: num_gto = 40
 
      INTEGER                                            :: i, ii, iw, k, l, ll, m, mb, mo, ngp, nn, &
                                                            nr, quadtype, relativistic, z
      INTEGER, DIMENSION(0:lmat)                         :: starti
      INTEGER, DIMENSION(0:lmat, 10)                     :: ncalc, ncore, nelem
      INTEGER, DIMENSION(:), POINTER                     :: econf
      LOGICAL                                            :: do_basopt, ecp_semi_local, monovalent
      REAL(kind=dp)                                      :: al, aval, cc, cval, ear, rk, xx, zeff
      REAL(kind=dp), DIMENSION(num_gto+2)                :: results
      TYPE(all_potential_type), POINTER                  :: all_potential
      TYPE(atom_basis_type), POINTER                     :: basis
      TYPE(atom_integrals), POINTER                      :: integrals
      TYPE(atom_orbitals), POINTER                       :: orbitals
      TYPE(atom_potential_type), POINTER                 :: potential
      TYPE(atom_type), POINTER                           :: atom
      TYPE(grid_atom_type), POINTER                      :: grid
      TYPE(gth_potential_type), POINTER                  :: gth_potential
      TYPE(sgp_potential_type), POINTER                  :: sgp_potential
 
      NULLIFY (atom)
      CALL create_atom_type(atom)
 
      CALL get_atomic_kind(atomic_kind, z=z)
      NULLIFY (all_potential, gth_potential)
      CALL get_qs_kind(qs_kind, zeff=zeff, &
                       all_potential=all_potential, &
                       gth_potential=gth_potential, &
                       sgp_potential=sgp_potential, &
                       monovalent=monovalent)
 
      IF (PRESENT(iunit)) THEN
         iw = iunit
      ELSE
         iw = -1
      END IF
 
      IF (PRESENT(allelectron)) THEN
         IF (allelectron) THEN
            NULLIFY (gth_potential)
            zeff = z
         END IF
      END IF
 
      do_basopt = .true.
      IF (PRESENT(optbasis)) THEN
         do_basopt = optbasis
      END IF
 
      cpassert(ngto <= num_gto)
 
      IF (ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)) THEN
         ! PP calculation are non-relativistic
         relativistic = do_nonrel_atom
      ELSE
         ! AE calculations use DKH2
         relativistic = do_dkh2_atom
      END IF
 
      atom%z = z
      CALL set_atom(atom, &
                    pp_calc=(ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)), &
                    method_type=do_rks_atom, &
                    relativistic=relativistic, &
                    coulomb_integral_type=do_numeric, &
                    exchange_integral_type=do_numeric)
 
      ALLOCATE (potential, basis, integrals)
 
      IF (PRESENT(confine)) THEN
         potential%confinement = confine
      ELSE
         IF (ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)) THEN
            potential%confinement = .true.
         ELSE
            potential%confinement = .false.
         END IF
      END IF
      potential%conf_type = barrier_conf
      potential%acon = 200._dp
      potential%rcon = 4.0_dp
      potential%scon = 8.0_dp
 
      IF (ASSOCIATED(gth_potential)) THEN
         potential%ppot_type = gth_pseudo
         CALL get_potential(gth_potential, zeff=zeff)
         CALL gth_potential_conversion(gth_potential, potential%gth_pot)
         CALL set_atom(atom, zcore=nint(zeff), potential=potential)
      ELSE IF (ASSOCIATED(sgp_potential)) THEN
         CALL get_potential(sgp_potential, ecp_semi_local=ecp_semi_local)
         IF (ALLOCATED(sgp_potential%upf_data)) THEN
            potential%ppot_type = upf_pseudo
            potential%upf_pot = sgp_potential%upf_data
         ELSE IF (ecp_semi_local) THEN
            potential%ppot_type = ecp_pseudo
            CALL ecp_potential_conversion(sgp_potential, potential%ecp_pot)
            potential%ecp_pot%symbol = ptable(z)%symbol
         ELSE
            potential%ppot_type = sgp_pseudo
            CALL sgp_potential_conversion(sgp_potential, potential%sgp_pot)
            potential%sgp_pot%symbol = ptable(z)%symbol
         END IF
         CALL get_potential(sgp_potential, zeff=zeff)
         CALL set_atom(atom, zcore=nint(zeff), potential=potential)
      ELSE
         potential%ppot_type = no_pseudo
         CALL set_atom(atom, zcore=z, potential=potential)
      END IF
 
      ! atomic grid
      NULLIFY (grid)
      ngp = 400
      quadtype = do_gapw_log
      CALL allocate_grid_atom(grid)
      CALL create_grid_atom(grid, ngp, 1, 1, 0, quadtype)
      grid%nr = ngp
      basis%grid => grid
 
      NULLIFY (basis%am, basis%cm, basis%as, basis%ns, basis%bf, basis%dbf, basis%ddbf)
 
      ! fill in the basis data structures
      basis%eps_eig = 1.e-12_dp
      basis%basis_type = gto_basis
      CALL clementi_geobas(z, cval, aval, basis%nbas, starti)
      basis%nprim = basis%nbas
      m = maxval(basis%nbas)
      ALLOCATE (basis%am(m, 0:lmat))
      basis%am = 0._dp
      DO l = 0, lmat
         DO i = 1, basis%nbas(l)
            ll = i - 1 + starti(l)
            basis%am(i, l) = aval*cval**(ll)
         END DO
      END DO
 
      basis%geometrical = .true.
      basis%aval = aval
      basis%cval = cval
      basis%start = starti
 
      ! initialize basis function on a radial grid
      nr = basis%grid%nr
      m = maxval(basis%nbas)
      ALLOCATE (basis%bf(nr, m, 0:lmat))
      ALLOCATE (basis%dbf(nr, m, 0:lmat))
      ALLOCATE (basis%ddbf(nr, m, 0:lmat))
      basis%bf = 0._dp
      basis%dbf = 0._dp
      basis%ddbf = 0._dp
      DO l = 0, lmat
         DO i = 1, basis%nbas(l)
            al = basis%am(i, l)
            DO k = 1, nr
               rk = basis%grid%rad(k)
               ear = exp(-al*basis%grid%rad(k)**2)
               basis%bf(k, i, l) = rk**l*ear
               basis%dbf(k, i, l) = (real(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear
               basis%ddbf(k, i, l) = (real(l*(l - 1), dp)*rk**(l - 2) - &
                                      2._dp*al*real(2*l + 1, dp)*rk**(l) + 4._dp*al**2*rk**(l + 2))*ear
            END DO
         END DO
      END DO
 
      CALL set_atom(atom, basis=basis)
 
      ! optimization defaults
      atom%optimization%damping = 0.2_dp
      atom%optimization%eps_scf = 1.e-6_dp
      atom%optimization%eps_diis = 100._dp
      atom%optimization%max_iter = 50
      atom%optimization%n_diis = 5
 
      nelem = 0
      ncore = 0
      ncalc = 0
      IF (monovalent) THEN
         ncalc(0, 1) = 1
         nelem(0, 1) = 1
      ELSE IF (ASSOCIATED(gth_potential)) THEN
         CALL get_potential(gth_potential, elec_conf=econf)
         CALL set_pseudo_state(econf, z, ncalc, ncore, nelem)
      ELSE IF (ASSOCIATED(sgp_potential)) THEN
         CALL get_potential(sgp_potential, elec_conf=econf)
         CALL set_pseudo_state(econf, z, ncalc, ncore, nelem)
      ELSE
         DO l = 0, min(lmat, ubound(ptable(z)%e_conv, 1))
            ll = 2*(2*l + 1)
            nn = ptable(z)%e_conv(l)
            ii = 0
            DO
               ii = ii + 1
               IF (nn <= ll) THEN
                  nelem(l, ii) = nn
                  EXIT
               ELSE
                  nelem(l, ii) = ll
                  nn = nn - ll
               END IF
            END DO
         END DO
         ncalc = nelem - ncore
      END IF
 
      IF (qs_kind%ghost .OR. qs_kind%floating) THEN
         nelem = 0
         ncore = 0
         ncalc = 0
      END IF
 
      ALLOCATE (atom%state)
 
      atom%state%core = 0._dp
      atom%state%core(0:lmat, 1:7) = real(ncore(0:lmat, 1:7), dp)
      atom%state%occ = 0._dp
      atom%state%occ(0:lmat, 1:7) = real(ncalc(0:lmat, 1:7), dp)
      atom%state%occupation = 0._dp
      atom%state%multiplicity = -1
      DO l = 0, lmat
         k = 0
         DO i = 1, 7
            IF (ncalc(l, i) > 0) THEN
               k = k + 1
               atom%state%occupation(l, k) = real(ncalc(l, i), dp)
            END IF
         END DO
      END DO
 
      atom%state%maxl_occ = get_maxl_occ(atom%state%occupation)
      atom%state%maxn_occ = get_maxn_occ(atom%state%occupation)
      atom%state%maxl_calc = atom%state%maxl_occ
      atom%state%maxn_calc = atom%state%maxn_occ
 
      ! calculate integrals
      ! general integrals
      CALL atom_int_setup(integrals, basis, potential=atom%potential, &
                          eri_coulomb=(atom%coulomb_integral_type == do_analytic), &
                          eri_exchange=(atom%exchange_integral_type == do_analytic))
      ! potential
      CALL atom_ppint_setup(integrals, basis, potential=atom%potential)
      ! relativistic correction terms
      NULLIFY (integrals%tzora, integrals%hdkh)
      CALL atom_relint_setup(integrals, basis, atom%relativistic, zcore=real(atom%zcore, dp))
      CALL set_atom(atom, integrals=integrals)
 
      NULLIFY (orbitals)
      mo = maxval(atom%state%maxn_calc)
      mb = maxval(atom%basis%nbas)
      CALL create_atom_orbs(orbitals, mb, mo)
      CALL set_atom(atom, orbitals=orbitals)
 
      CALL calculate_atom(atom, iw)
 
      IF (do_basopt) THEN
         CALL atom_fit_density(atom, ngto, 0, iw, results=results)
         xx = results(1)
         cc = results(2)
         DO i = 1, ngto
            density(i, 1) = xx*cc**i
            density(i, 2) = results(2 + i)
         END DO
      ELSE
         CALL atom_fit_density(atom, ngto, 0, iw, agto=density(:, 1), results=results)
         density(1:ngto, 2) = results(1:ngto)
      END IF
 
      ! clean up
      CALL atom_int_release(integrals)
      CALL atom_ppint_release(integrals)
      CALL atom_relint_release(integrals)
      CALL release_atom_basis(basis)
      CALL release_atom_potential(potential)
      CALL release_atom_type(atom)
 
      DEALLOCATE (potential, basis, integrals)
 
   END SUBROUTINE calculate_atomic_density
 
! **************************************************************************************************
!> \brief ...
!> \param atomic_kind ...
!> \param qs_kind ...
!> \param rel_control ...
!> \param rtmat ...
! **************************************************************************************************
   SUBROUTINE calculate_atomic_relkin(atomic_kind, qs_kind, rel_control, rtmat)
      TYPE(atomic_kind_type), INTENT(IN)                 :: atomic_kind
      TYPE(qs_kind_type), INTENT(IN)                     :: qs_kind
      TYPE(rel_control_type), POINTER                    :: rel_control
      REAL(kind=dp), DIMENSION(:, :), POINTER            :: rtmat
 
      INTEGER                                            :: i, ii, ipgf, j, k, k1, k2, l, ll, m, n, &
                                                            ngp, nj, nn, nr, ns, nset, nsgf, &
                                                            quadtype, relativistic, z
      INTEGER, DIMENSION(0:lmat, 10)                     :: ncalc, ncore, nelem
      INTEGER, DIMENSION(0:lmat, 100)                    :: set_index, shell_index
      INTEGER, DIMENSION(:), POINTER                     :: lmax, lmin, npgf, nshell
      INTEGER, DIMENSION(:, :), POINTER                  :: first_sgf, last_sgf, ls
      REAL(kind=dp)                                      :: al, alpha, ear, prefac, rk, zeff
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: omat
      REAL(kind=dp), DIMENSION(:, :), POINTER            :: zet
      REAL(kind=dp), DIMENSION(:, :, :), POINTER         :: gcc
      TYPE(all_potential_type), POINTER                  :: all_potential
      TYPE(atom_basis_type), POINTER                     :: basis
      TYPE(atom_integrals), POINTER                      :: integrals
      TYPE(atom_potential_type), POINTER                 :: potential
      TYPE(atom_type), POINTER                           :: atom
      TYPE(grid_atom_type), POINTER                      :: grid
      TYPE(gto_basis_set_type), POINTER                  :: orb_basis_set
 
      IF (rel_control%rel_method == rel_none) RETURN
 
      NULLIFY (all_potential, orb_basis_set)
      CALL get_qs_kind(qs_kind, basis_set=orb_basis_set, all_potential=all_potential)
 
      cpassert(ASSOCIATED(orb_basis_set))
 
      IF (ASSOCIATED(all_potential)) THEN
         ! only all electron atoms will get the relativistic correction
 
         CALL get_atomic_kind(atomic_kind, z=z)
         CALL get_qs_kind(qs_kind, zeff=zeff)
         NULLIFY (atom)
         CALL create_atom_type(atom)
         NULLIFY (atom%xc_section)
         NULLIFY (atom%orbitals)
         atom%z = z
         alpha = sqrt(all_potential%alpha_core_charge)
 
         ! set the method flag
         SELECT CASE (rel_control%rel_method)
         CASE DEFAULT
            cpabort("Unknown relativistic method for calculate_atomic_relkin")
         CASE (rel_dkh)
            SELECT CASE (rel_control%rel_DKH_order)
            CASE DEFAULT
               cpabort("Unknown DKH order for calculate_atomic_relkin")
            CASE (0)
               relativistic = do_dkh0_atom
            CASE (1)
               relativistic = do_dkh1_atom
            CASE (2)
               relativistic = do_dkh2_atom
            CASE (3)
               relativistic = do_dkh3_atom
            END SELECT
         CASE (rel_zora)
            SELECT CASE (rel_control%rel_zora_type)
            CASE DEFAULT
               cpabort("Unknown ZORA type for calculate_atomic_relkin")
            CASE (rel_zora_full)
               cpabort("ZORA full not yet implemented in calculate_atomic_relkin")
            CASE (rel_zora_mp)
               relativistic = do_zoramp_atom
            CASE (rel_sczora_mp)
               relativistic = do_sczoramp_atom
            END SELECT
         END SELECT
 
         CALL set_atom(atom, &
                       pp_calc=.false., &
                       method_type=do_rks_atom, &
                       relativistic=relativistic, &
                       coulomb_integral_type=do_numeric, &
                       exchange_integral_type=do_numeric)
 
         ALLOCATE (potential, basis, integrals)
 
         potential%ppot_type = no_pseudo
         CALL set_atom(atom, zcore=z, potential=potential)
 
         CALL get_gto_basis_set(orb_basis_set, &
                                nset=nset, nshell=nshell, npgf=npgf, lmin=lmin, lmax=lmax, l=ls, nsgf=nsgf, zet=zet, gcc=gcc, &
                                first_sgf=first_sgf, last_sgf=last_sgf)
 
         NULLIFY (grid)
         ngp = 400
         quadtype = do_gapw_log
         CALL allocate_grid_atom(grid)
         CALL create_grid_atom(grid, ngp, 1, 1, 0, quadtype)
         grid%nr = ngp
         basis%grid => grid
 
         NULLIFY (basis%am, basis%cm, basis%as, basis%ns, basis%bf, basis%dbf, basis%ddbf)
         basis%basis_type = cgto_basis
         basis%eps_eig = 1.e-12_dp
 
         ! fill in the basis data structures
         set_index = 0
         shell_index = 0
         basis%nprim = 0
         basis%nbas = 0
         DO i = 1, nset
            DO j = lmin(i), min(lmax(i), lmat)
               basis%nprim(j) = basis%nprim(j) + npgf(i)
            END DO
            DO j = 1, nshell(i)
               l = ls(j, i)
               IF (l <= lmat) THEN
                  basis%nbas(l) = basis%nbas(l) + 1
                  k = basis%nbas(l)
                  cpassert(k <= 100)
                  set_index(l, k) = i
                  shell_index(l, k) = j
               END IF
            END DO
         END DO
 
         nj = maxval(basis%nprim)
         ns = maxval(basis%nbas)
         ALLOCATE (basis%am(nj, 0:lmat))
         basis%am = 0._dp
         ALLOCATE (basis%cm(nj, ns, 0:lmat))
         basis%cm = 0._dp
         DO j = 0, lmat
            nj = 0
            ns = 0
            DO i = 1, nset
               IF (j >= lmin(i) .AND. j <= lmax(i)) THEN
                  DO ipgf = 1, npgf(i)
                     basis%am(nj + ipgf, j) = zet(ipgf, i)
                  END DO
                  DO ii = 1, nshell(i)
                     IF (ls(ii, i) == j) THEN
                        ns = ns + 1
                        DO ipgf = 1, npgf(i)
                           basis%cm(nj + ipgf, ns, j) = gcc(ipgf, ii, i)
                        END DO
                     END IF
                  END DO
                  nj = nj + npgf(i)
               END IF
            END DO
         END DO
 
         ! Normalization as used in the atomic code
         ! We have to undo the Quickstep normalization
         DO j = 0, lmat
            prefac = 2.0_dp*sqrt(pi/dfac(2*j + 1))
            DO ipgf = 1, basis%nprim(j)
               DO ii = 1, basis%nbas(j)
                  basis%cm(ipgf, ii, j) = prefac*basis%cm(ipgf, ii, j)
               END DO
            END DO
         END DO
 
         ! initialize basis function on a radial grid
         nr = basis%grid%nr
         m = maxval(basis%nbas)
         ALLOCATE (basis%bf(nr, m, 0:lmat))
         ALLOCATE (basis%dbf(nr, m, 0:lmat))
         ALLOCATE (basis%ddbf(nr, m, 0:lmat))
 
         basis%bf = 0._dp
         basis%dbf = 0._dp
         basis%ddbf = 0._dp
         DO l = 0, lmat
            DO i = 1, basis%nprim(l)
               al = basis%am(i, l)
               DO k = 1, nr
                  rk = basis%grid%rad(k)
                  ear = exp(-al*basis%grid%rad(k)**2)
                  DO j = 1, basis%nbas(l)
                     basis%bf(k, j, l) = basis%bf(k, j, l) + rk**l*ear*basis%cm(i, j, l)
                     basis%dbf(k, j, l) = basis%dbf(k, j, l) &
                                          + (real(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear*basis%cm(i, j, l)
                     basis%ddbf(k, j, l) = basis%ddbf(k, j, l) + &
                                           (real(l*(l - 1), dp)*rk**(l - 2) - 2._dp*al*real(2*l + 1, dp)* &
                                            rk**(l) + 4._dp*al**2*rk**(l + 2))*ear*basis%cm(i, j, l)
                  END DO
               END DO
            END DO
         END DO
 
         CALL set_atom(atom, basis=basis)
 
         ! optimization defaults
         atom%optimization%damping = 0.2_dp
         atom%optimization%eps_scf = 1.e-6_dp
         atom%optimization%eps_diis = 100._dp
         atom%optimization%max_iter = 50
         atom%optimization%n_diis = 5
 
         ! electronic state
         nelem = 0
         ncore = 0
         ncalc = 0
         DO l = 0, min(lmat, ubound(ptable(z)%e_conv, 1))
            ll = 2*(2*l + 1)
            nn = ptable(z)%e_conv(l)
            ii = 0
            DO
               ii = ii + 1
               IF (nn <= ll) THEN
                  nelem(l, ii) = nn
                  EXIT
               ELSE
                  nelem(l, ii) = ll
                  nn = nn - ll
               END IF
            END DO
         END DO
         ncalc = nelem - ncore
 
         IF (qs_kind%ghost .OR. qs_kind%floating) THEN
            nelem = 0
            ncore = 0
            ncalc = 0
         END IF
 
         ALLOCATE (atom%state)
 
         atom%state%core = 0._dp
         atom%state%core(0:lmat, 1:7) = real(ncore(0:lmat, 1:7), dp)
         atom%state%occ = 0._dp
         atom%state%occ(0:lmat, 1:7) = real(ncalc(0:lmat, 1:7), dp)
         atom%state%occupation = 0._dp
         atom%state%multiplicity = -1
         DO l = 0, lmat
            k = 0
            DO i = 1, 7
               IF (ncalc(l, i) > 0) THEN
                  k = k + 1
                  atom%state%occupation(l, k) = real(ncalc(l, i), dp)
               END IF
            END DO
         END DO
 
         atom%state%maxl_occ = get_maxl_occ(atom%state%occupation)
         atom%state%maxn_occ = get_maxn_occ(atom%state%occupation)
         atom%state%maxl_calc = atom%state%maxl_occ
         atom%state%maxn_calc = atom%state%maxn_occ
 
         ! calculate integrals
         ! general integrals
         CALL atom_int_setup(integrals, basis)
         ! potential
         CALL atom_ppint_setup(integrals, basis, potential=atom%potential)
         ! relativistic correction terms
         NULLIFY (integrals%tzora, integrals%hdkh)
         CALL atom_relint_setup(integrals, basis, atom%relativistic, zcore=real(atom%zcore, dp), &
                                alpha=alpha)
         CALL set_atom(atom, integrals=integrals)
 
         ! for DKH we need erfc integrals to correct non-relativistic
         integrals%core = 0.0_dp
         DO l = 0, lmat
            n = integrals%n(l)
            m = basis%nprim(l)
            ALLOCATE (omat(m, m))
 
            CALL sg_erfc(omat(1:m, 1:m), l, alpha, basis%am(1:m, l), basis%am(1:m, l))
            integrals%core(1:n, 1:n, l) = matmul(transpose(basis%cm(1:m, 1:n, l)), &
                                                 matmul(omat(1:m, 1:m), basis%cm(1:m, 1:n, l)))
 
            DEALLOCATE (omat)
         END DO
 
         ! recover relativistic kinetic matrix in CP2K/GPW order and normalization
         IF (ASSOCIATED(rtmat)) THEN
            DEALLOCATE (rtmat)
         END IF
         ALLOCATE (rtmat(nsgf, nsgf))
         rtmat = 0._dp
         DO l = 0, lmat
            ll = 2*l
            DO k1 = 1, basis%nbas(l)
               DO k2 = 1, basis%nbas(l)
                  i = first_sgf(shell_index(l, k1), set_index(l, k1))
                  j = first_sgf(shell_index(l, k2), set_index(l, k2))
                  SELECT CASE (atom%relativistic)
                  CASE DEFAULT
                     cpabort("Unknown relativistic type for calculate_atomic_relkin")
                  CASE (do_zoramp_atom, do_sczoramp_atom)
                     DO m = 0, ll
                        rtmat(i + m, j + m) = integrals%tzora(k1, k2, l)
                     END DO
                  CASE (do_dkh0_atom, do_dkh1_atom, do_dkh2_atom, do_dkh3_atom)
                     DO m = 0, ll
                        rtmat(i + m, j + m) = integrals%hdkh(k1, k2, l) - integrals%kin(k1, k2, l) + &
                                              atom%zcore*integrals%core(k1, k2, l)
                     END DO
                  END SELECT
               END DO
            END DO
         END DO
         DO k1 = 1, nsgf
            DO k2 = k1, nsgf
               rtmat(k1, k2) = 0.5_dp*(rtmat(k1, k2) + rtmat(k2, k1))
               rtmat(k2, k1) = rtmat(k1, k2)
            END DO
         END DO
 
         ! clean up
         CALL atom_int_release(integrals)
         CALL atom_ppint_release(integrals)
         CALL atom_relint_release(integrals)
         CALL release_atom_basis(basis)
         CALL release_atom_potential(potential)
         CALL release_atom_type(atom)
 
         DEALLOCATE (potential, basis, integrals)
 
      ELSE
 
         IF (ASSOCIATED(rtmat)) THEN
            DEALLOCATE (rtmat)
         END IF
         NULLIFY (rtmat)
 
      END IF
 
   END SUBROUTINE calculate_atomic_relkin
 
! **************************************************************************************************
!> \brief ...
!> \param gth_potential ...
!> \param gth_atompot ...
! **************************************************************************************************
   SUBROUTINE gth_potential_conversion(gth_potential, gth_atompot)
      TYPE(gth_potential_type), POINTER                  :: gth_potential
      TYPE(atom_gthpot_type)                             :: gth_atompot
 
      INTEGER                                            :: i, j, l, lm, n, ne, nexp_lpot, nexp_lsd, &
                                                            nexp_nlcc
      INTEGER, DIMENSION(:), POINTER                     :: nct_lpot, nct_lsd, nct_nlcc, nppnl, &
                                                            ppeconf
      LOGICAL                                            :: lpot_present, lsd_present, nlcc_present, &
                                                            soc_present
      REAL(kind=dp)                                      :: ac, zeff
      REAL(kind=dp), DIMENSION(:), POINTER               :: alpha_lpot, alpha_lsd, alpha_nlcc, ap, ce
      REAL(kind=dp), DIMENSION(:, :), POINTER            :: cval_lpot, cval_lsd, cval_nlcc
      REAL(kind=dp), DIMENSION(:, :, :), POINTER         :: hp, kp
 
      CALL get_potential(gth_potential, &
                         zeff=zeff, &
                         elec_conf=ppeconf, &
                         alpha_core_charge=ac, &
                         nexp_ppl=ne, &
                         cexp_ppl=ce, &
                         lppnl=lm, &
                         nprj_ppnl=nppnl, &
                         alpha_ppnl=ap, &
                         kprj_ppnl=kp, &
                         hprj_ppnl=hp)
 
      gth_atompot%zion = zeff
      gth_atompot%rc = sqrt(0.5_dp/ac)
      gth_atompot%ncl = ne
      gth_atompot%cl(:) = 0._dp
      IF (ac > 0._dp) THEN
         DO i = 1, ne
            gth_atompot%cl(i) = ce(i)/(2._dp*ac)**(i - 1)
         END DO
      END IF
      !extended type
      gth_atompot%lpotextended = .false.
      gth_atompot%lsdpot = .false.
      gth_atompot%nlcc = .false.
      gth_atompot%nexp_lpot = 0
      gth_atompot%nexp_lsd = 0
      gth_atompot%nexp_nlcc = 0
      CALL get_potential(gth_potential, &
                         lpot_present=lpot_present, &
                         lsd_present=lsd_present, &
                         nlcc_present=nlcc_present)
      IF (lpot_present) THEN
         CALL get_potential(gth_potential, &
                            nexp_lpot=nexp_lpot, &
                            alpha_lpot=alpha_lpot, &
                            nct_lpot=nct_lpot, &
                            cval_lpot=cval_lpot)
         gth_atompot%lpotextended = .true.
         gth_atompot%nexp_lpot = nexp_lpot
         gth_atompot%alpha_lpot(1:nexp_lpot) = sqrt(0.5_dp/alpha_lpot(1:nexp_lpot))
         gth_atompot%nct_lpot(1:nexp_lpot) = nct_lpot(1:nexp_lpot)
         DO j = 1, nexp_lpot
            ac = alpha_lpot(j)
            DO i = 1, 4
               gth_atompot%cval_lpot(i, j) = cval_lpot(i, j)/(2._dp*ac)**(i - 1)
            END DO
         END DO
      END IF
      IF (lsd_present) THEN
         CALL get_potential(gth_potential, &
                            nexp_lsd=nexp_lsd, &
                            alpha_lsd=alpha_lsd, &
                            nct_lsd=nct_lsd, &
                            cval_lsd=cval_lsd)
         gth_atompot%lsdpot = .true.
         gth_atompot%nexp_lsd = nexp_lsd
         gth_atompot%alpha_lsd(1:nexp_lsd) = sqrt(0.5_dp/alpha_lsd(1:nexp_lsd))
         gth_atompot%nct_lsd(1:nexp_lsd) = nct_lsd(1:nexp_lsd)
         DO j = 1, nexp_lpot
            ac = alpha_lsd(j)
            DO i = 1, 4
               gth_atompot%cval_lsd(i, j) = cval_lsd(i, j)/(2._dp*ac)**(i - 1)
            END DO
         END DO
      END IF
 
      ! nonlocal part
      gth_atompot%nl(:) = 0
      gth_atompot%rcnl(:) = 0._dp
      gth_atompot%hnl(:, :, :) = 0._dp
      DO l = 0, lm
         n = nppnl(l)
         gth_atompot%nl(l) = n
         gth_atompot%rcnl(l) = sqrt(0.5_dp/ap(l))
         gth_atompot%hnl(1:n, 1:n, l) = hp(1:n, 1:n, l)
      END DO
 
      ! SOC
      CALL get_potential(gth_potential, soc_present=soc_present)
      gth_atompot%soc = soc_present
      gth_atompot%knl = 0.0_dp
      IF (soc_present) THEN
         DO l = 1, lm
            n = nppnl(l)
            gth_atompot%knl(1:n, 1:n, l) = kp(1:n, 1:n, l)
         END DO
      END IF
 
      IF (nlcc_present) THEN
         CALL get_potential(gth_potential, &
                            nexp_nlcc=nexp_nlcc, &
                            alpha_nlcc=alpha_nlcc, &
                            nct_nlcc=nct_nlcc, &
                            cval_nlcc=cval_nlcc)
         gth_atompot%nlcc = .true.
         gth_atompot%nexp_nlcc = nexp_nlcc
         gth_atompot%alpha_nlcc(1:nexp_nlcc) = alpha_nlcc(1:nexp_nlcc)
         gth_atompot%nct_nlcc(1:nexp_nlcc) = nct_nlcc(1:nexp_nlcc)
         gth_atompot%cval_nlcc(1:4, 1:nexp_nlcc) = cval_nlcc(1:4, 1:nexp_nlcc)
      END IF
 
   END SUBROUTINE gth_potential_conversion
 
! **************************************************************************************************
!> \brief ...
!> \param sgp_potential ...
!> \param sgp_atompot ...
! **************************************************************************************************
   SUBROUTINE sgp_potential_conversion(sgp_potential, sgp_atompot)
      TYPE(sgp_potential_type), POINTER                  :: sgp_potential
      TYPE(atom_sgppot_type)                             :: sgp_atompot
 
      INTEGER                                            :: lm, n
      INTEGER, DIMENSION(:), POINTER                     :: ppeconf
      LOGICAL                                            :: nlcc_present
      REAL(kind=dp)                                      :: ac, zeff
      REAL(kind=dp), DIMENSION(:), POINTER               :: ap, ce
      REAL(kind=dp), DIMENSION(:, :), POINTER            :: hhp
      REAL(kind=dp), DIMENSION(:, :, :), POINTER         :: ccp
 
      CALL get_potential(sgp_potential, &
                         name=sgp_atompot%pname, &
                         zeff=zeff, &
                         elec_conf=ppeconf, &
                         core_charge_radius=ac)
      sgp_atompot%zion = zeff
      sgp_atompot%ac_local = ac
      sgp_atompot%econf(0:3) = ppeconf(0:3)
      CALL get_potential(sgp_potential, lmax=lm, &
                         is_nonlocal=sgp_atompot%is_nonlocal, &
                         n_nonlocal=n, a_nonlocal=ap, h_nonlocal=hhp, c_nonlocal=ccp)
      ! nonlocal
      sgp_atompot%has_nonlocal = any(sgp_atompot%is_nonlocal)
      sgp_atompot%lmax = lm
      IF (sgp_atompot%has_nonlocal) THEN
         cpassert(n <= SIZE(sgp_atompot%a_nonlocal))
         sgp_atompot%n_nonlocal = n
         sgp_atompot%a_nonlocal(1:n) = ap(1:n)
         sgp_atompot%h_nonlocal(1:n, 0:lm) = hhp(1:n, 0:lm)
         ! The atomic solver uses unit-normalized spherical harmonics and
         ! radial integrals; undo the SGP storage factor 1/sqrt(4*pi).
         sgp_atompot%c_nonlocal(1:n, 1:n, 0:lm) = sqrt(fourpi)*ccp(1:n, 1:n, 0:lm)
      END IF
      ! local
      CALL get_potential(sgp_potential, n_local=n, a_local=ap, c_local=ce)
      cpassert(n <= SIZE(sgp_atompot%a_local))
      sgp_atompot%n_local = n
      sgp_atompot%a_local(1:n) = ap(1:n)
      sgp_atompot%c_local(1:n) = ce(1:n)
      ! NLCC
      CALL get_potential(sgp_potential, has_nlcc=nlcc_present, &
                         n_nlcc=n, a_nlcc=ap, c_nlcc=ce)
      IF (nlcc_present) THEN
         sgp_atompot%has_nlcc = .true.
         cpassert(n <= SIZE(sgp_atompot%a_nlcc))
         sgp_atompot%n_nlcc = n
         sgp_atompot%a_nlcc(1:n) = ap(1:n)
         sgp_atompot%c_nlcc(1:n) = ce(1:n)
      ELSE
         sgp_atompot%has_nlcc = .false.
      END IF
 
   END SUBROUTINE sgp_potential_conversion
 
! **************************************************************************************************
!> \brief ...
!> \param sgp_potential ...
!> \param ecp_atompot ...
! **************************************************************************************************
   SUBROUTINE ecp_potential_conversion(sgp_potential, ecp_atompot)
      TYPE(sgp_potential_type), POINTER                  :: sgp_potential
      TYPE(atom_ecppot_type)                             :: ecp_atompot
 
      INTEGER, DIMENSION(:), POINTER                     :: ppeconf
      LOGICAL                                            :: ecp_local, ecp_semi_local
      REAL(kind=dp)                                      :: zeff
 
      CALL get_potential(sgp_potential, ecp_local=ecp_local, ecp_semi_local=ecp_semi_local)
      cpassert(ecp_semi_local .AND. ecp_local)
      CALL get_potential(sgp_potential, &
                         name=ecp_atompot%pname, &
                         zeff=zeff, &
                         elec_conf=ppeconf)
      ecp_atompot%zion = zeff
      ecp_atompot%econf(0:3) = ppeconf(0:3)
      CALL get_potential(sgp_potential, sl_lmax=ecp_atompot%lmax)
      ! local
      CALL get_potential(sgp_potential, nloc=ecp_atompot%nloc, nrloc=ecp_atompot%nrloc, &
                         aloc=ecp_atompot%aloc, bloc=ecp_atompot%bloc)
      ! nonlocal
      CALL get_potential(sgp_potential, npot=ecp_atompot%npot, nrpot=ecp_atompot%nrpot, &
                         apot=ecp_atompot%apot, bpot=ecp_atompot%bpot)
 
   END SUBROUTINE ecp_potential_conversion
! **************************************************************************************************
 
END MODULE atom_kind_orbitals