!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2020 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!

! **************************************************************************************************
!> \brief   Calculate the atomic operator matrices
!> \author  jgh
!> \date    03.03.2008
!> \version 1.0
!>
! **************************************************************************************************
MODULE atom_operators
   USE ai_onecenter,                    ONLY: &
        sg_coulomb, sg_erf, sg_erfc, sg_exchange, sg_gpot, sg_kinetic, sg_kinnuc, sg_nuclear, &
        sg_overlap, sg_proj_ol, sto_kinetic, sto_nuclear, sto_overlap
   USE atom_types,                      ONLY: &
        atom_basis_gridrep, atom_basis_type, atom_compare_grids, atom_integrals, &
        atom_potential_type, atom_state, cgto_basis, ecp_pseudo, gth_pseudo, gto_basis, lmat, &
        no_pseudo, num_basis, release_atom_basis, sgp_pseudo, sto_basis, upf_pseudo
   USE atom_utils,                      ONLY: &
        atom_solve, contract2, contract2add, contract4, coulomb_potential_numeric, integrate_grid, &
        numpot_matrix, slater_density, wigner_slater_functional
   USE dkh_main,                        ONLY: dkh_atom_transformation
   USE input_constants,                 ONLY: &
        barrier_conf, do_dkh0_atom, do_dkh1_atom, do_dkh2_atom, do_dkh3_atom, do_nonrel_atom, &
        do_sczoramp_atom, do_zoramp_atom, poly_conf
   USE kinds,                           ONLY: dp
   USE lapack,                          ONLY: lapack_sgesv
   USE mathconstants,                   ONLY: gamma1,&
                                              sqrt2
   USE mathlib,                         ONLY: jacobi
   USE periodic_table,                  ONLY: ptable
   USE physcon,                         ONLY: c_light_au
   USE qs_grid_atom,                    ONLY: grid_atom_type
#include "./base/base_uses.f90"

   IMPLICIT NONE

   PRIVATE

   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_operators'

   PUBLIC :: atom_int_setup, atom_ppint_setup, atom_int_release, atom_ppint_release
   PUBLIC :: atom_relint_setup, atom_relint_release, atom_basis_projection_overlap
   PUBLIC :: calculate_model_potential

CONTAINS

! **************************************************************************************************
!> \brief Set up atomic integrals.
!> \param integrals     atomic integrals
!> \param basis         atomic basis set
!> \param potential     pseudo-potential
!> \param eri_coulomb   setup one-centre Coulomb Electron Repulsion Integrals
!> \param eri_exchange  setup one-centre exchange Electron Repulsion Integrals
!> \param all_nu        compute integrals for all even integer parameters [0 .. 2*l]
!>                      REDUNDANT, AS THIS SUBROUTINE IS NEVER INVOKED WITH all_nu = .TRUE.
! **************************************************************************************************
   SUBROUTINE atom_int_setup(integrals, basis, potential, &
                             eri_coulomb, eri_exchange, all_nu)
      TYPE(atom_integrals), INTENT(INOUT)                :: integrals
      TYPE(atom_basis_type), INTENT(INOUT)               :: basis
      TYPE(atom_potential_type), INTENT(IN), OPTIONAL    :: potential
      LOGICAL, INTENT(IN), OPTIONAL                      :: eri_coulomb, eri_exchange, all_nu

      CHARACTER(len=*), PARAMETER                        :: routineN = 'atom_int_setup'

      INTEGER                                            :: handle, i, ii, info, ipiv(1000), l, l1, &
                                                            l2, ll, lwork, m, m1, m2, mm1, mm2, n, &
                                                            n1, n2, nn1, nn2, nu, nx
      REAL(KIND=dp)                                      :: om, rc, ron, sc, x
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: cpot, w, work
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: omat, vmat
      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: eri

      CALL timeset(routineN, handle)

      IF (integrals%status == 0) THEN
         n = MAXVAL(basis%nbas)
         integrals%n = basis%nbas

         IF (PRESENT(eri_coulomb)) THEN
            integrals%eri_coulomb = eri_coulomb
         ELSE
            integrals%eri_coulomb = .FALSE.
         END IF
         IF (PRESENT(eri_exchange)) THEN
            integrals%eri_exchange = eri_exchange
         ELSE
            integrals%eri_exchange = .FALSE.
         END IF
         IF (PRESENT(all_nu)) THEN
            integrals%all_nu = all_nu
         ELSE
            integrals%all_nu = .FALSE.
         END IF

         NULLIFY (integrals%ovlp, integrals%kin, integrals%core, integrals%conf)
         DO ll = 1, SIZE(integrals%ceri)
            NULLIFY (integrals%ceri(ll)%int, integrals%eeri(ll)%int)
         END DO

         ALLOCATE (integrals%ovlp(n, n, 0:lmat))
         integrals%ovlp = 0._dp

         ALLOCATE (integrals%kin(n, n, 0:lmat))
         integrals%kin = 0._dp

         integrals%status = 1

         IF (PRESENT(potential)) THEN
            IF (potential%confinement) THEN
               ALLOCATE (integrals%conf(n, n, 0:lmat))
               integrals%conf = 0._dp
               m = basis%grid%nr
               ALLOCATE (cpot(1:m))
               IF (potential%conf_type == poly_conf) THEN
                  rc = potential%rcon
                  sc = potential%scon
                  cpot(1:m) = (basis%grid%rad(1:m)/rc)**sc
               ELSEIF (potential%conf_type == barrier_conf) THEN
                  om = potential%rcon
                  ron = potential%scon
                  rc = ron + om
                  DO i = 1, m
                     IF (basis%grid%rad(i) < ron) THEN
                        cpot(i) = 0.0_dp
                     ELSEIF (basis%grid%rad(i) < rc) THEN
                        x = (basis%grid%rad(i) - ron)/om
                        x = 1._dp - x
                        cpot(i) = -6._dp*x**5 + 15._dp*x**4 - 10._dp*x**3 + 1._dp
                        x = (rc - basis%grid%rad(i))**2/om/(basis%grid%rad(i) - ron)
                        cpot(i) = cpot(i)*x
                     ELSE
                        cpot(i) = 1.0_dp
                     END IF
                  END DO
               ELSE
                  CPABORT("")
               END IF
               CALL numpot_matrix(integrals%conf, cpot, basis, 0)
               DEALLOCATE (cpot)
            END IF
         END IF

         SELECT CASE (basis%basis_type)
         CASE DEFAULT
            CPABORT("")
         CASE (GTO_BASIS)
            DO l = 0, lmat
               n = integrals%n(l)
               CALL sg_overlap(integrals%ovlp(1:n, 1:n, l), l, basis%am(1:n, l), basis%am(1:n, l))
               CALL sg_kinetic(integrals%kin(1:n, 1:n, l), l, basis%am(1:n, l), basis%am(1:n, l))
            END DO
            IF (integrals%eri_coulomb) THEN
               ll = 0
               DO l1 = 0, lmat
                  n1 = integrals%n(l1)
                  nn1 = (n1*(n1 + 1))/2
                  DO l2 = 0, l1
                     n2 = integrals%n(l2)
                     nn2 = (n2*(n2 + 1))/2
                     IF (integrals%all_nu) THEN
                        nx = MIN(2*l1, 2*l2)
                     ELSE
                        nx = 0
                     END IF
                     DO nu = 0, nx, 2
                        ll = ll + 1
                        CPASSERT(ll <= SIZE(integrals%ceri))
                        ALLOCATE (integrals%ceri(ll)%int(nn1, nn2))
                        integrals%ceri(ll)%int = 0._dp
                        eri => integrals%ceri(ll)%int
                        CALL sg_coulomb(eri, nu, basis%am(1:n1, l1), l1, basis%am(1:n2, l2), l2)
                     END DO
                  END DO
               END DO
            END IF
            IF (integrals%eri_exchange) THEN
               ll = 0
               DO l1 = 0, lmat
                  n1 = integrals%n(l1)
                  nn1 = (n1*(n1 + 1))/2
                  DO l2 = 0, l1
                     n2 = integrals%n(l2)
                     nn2 = (n2*(n2 + 1))/2
                     DO nu = ABS(l1 - l2), l1 + l2, 2
                        ll = ll + 1
                        CPASSERT(ll <= SIZE(integrals%eeri))
                        ALLOCATE (integrals%eeri(ll)%int(nn1, nn2))
                        integrals%eeri(ll)%int = 0._dp
                        eri => integrals%eeri(ll)%int
                        CALL sg_exchange(eri, nu, basis%am(1:n1, l1), l1, basis%am(1:n2, l2), l2)
                     END DO
                  END DO
               END DO
            END IF
         CASE (CGTO_BASIS)
            DO l = 0, lmat
               n = integrals%n(l)
               m = basis%nprim(l)
               IF (n > 0 .AND. m > 0) THEN
                  ALLOCATE (omat(m, m))

                  CALL sg_overlap(omat(1:m, 1:m), l, basis%am(1:m, l), basis%am(1:m, l))
                  CALL contract2(integrals%ovlp(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l))
                  CALL sg_kinetic(omat(1:m, 1:m), l, basis%am(1:m, l), basis%am(1:m, l))
                  CALL contract2(integrals%kin(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l))
                  DEALLOCATE (omat)
               END IF
            END DO
            IF (integrals%eri_coulomb) THEN
               ll = 0
               DO l1 = 0, lmat
                  n1 = integrals%n(l1)
                  nn1 = (n1*(n1 + 1))/2
                  m1 = basis%nprim(l1)
                  mm1 = (m1*(m1 + 1))/2
                  DO l2 = 0, l1
                     n2 = integrals%n(l2)
                     nn2 = (n2*(n2 + 1))/2
                     m2 = basis%nprim(l2)
                     mm2 = (m2*(m2 + 1))/2
                     IF (integrals%all_nu) THEN
                        nx = MIN(2*l1, 2*l2)
                     ELSE
                        nx = 0
                     END IF
                     DO nu = 0, nx, 2
                        ll = ll + 1
                        CPASSERT(ll <= SIZE(integrals%ceri))
                        ALLOCATE (integrals%ceri(ll)%int(nn1, nn2))
                        integrals%ceri(ll)%int = 0._dp
                        ALLOCATE (omat(mm1, mm2))

                        eri => integrals%ceri(ll)%int
                        CALL sg_coulomb(omat, nu, basis%am(1:m1, l1), l1, basis%am(1:m2, l2), l2)
                        CALL contract4(eri, omat, basis%cm(1:m1, 1:n1, l1), basis%cm(1:m2, 1:n2, l2))

                        DEALLOCATE (omat)
                     END DO
                  END DO
               END DO
            END IF
            IF (integrals%eri_exchange) THEN
               ll = 0
               DO l1 = 0, lmat
                  n1 = integrals%n(l1)
                  nn1 = (n1*(n1 + 1))/2
                  m1 = basis%nprim(l1)
                  mm1 = (m1*(m1 + 1))/2
                  DO l2 = 0, l1
                     n2 = integrals%n(l2)
                     nn2 = (n2*(n2 + 1))/2
                     m2 = basis%nprim(l2)
                     mm2 = (m2*(m2 + 1))/2
                     DO nu = ABS(l1 - l2), l1 + l2, 2
                        ll = ll + 1
                        CPASSERT(ll <= SIZE(integrals%eeri))
                        ALLOCATE (integrals%eeri(ll)%int(nn1, nn2))
                        integrals%eeri(ll)%int = 0._dp
                        ALLOCATE (omat(mm1, mm2))

                        eri => integrals%eeri(ll)%int
                        CALL sg_exchange(omat, nu, basis%am(1:m1, l1), l1, basis%am(1:m2, l2), l2)
                        CALL contract4(eri, omat, basis%cm(1:m1, 1:n1, l1), basis%cm(1:m2, 1:n2, l2))

                        DEALLOCATE (omat)
                     END DO
                  END DO
               END DO
            END IF
         CASE (STO_BASIS)
            DO l = 0, lmat
               n = integrals%n(l)
               CALL sto_overlap(integrals%ovlp(1:n, 1:n, l), basis%ns(1:n, l), basis%as(1:n, l), &
                                basis%ns(1:n, l), basis%as(1:n, l))
               CALL sto_kinetic(integrals%kin(1:n, 1:n, l), l, basis%ns(1:n, l), basis%as(1:n, l), &
                                basis%ns(1:n, l), basis%as(1:n, l))
            END DO
            CPASSERT(.NOT. integrals%eri_coulomb)
            CPASSERT(.NOT. integrals%eri_exchange)
         CASE (NUM_BASIS)
            CPABORT("")
         END SELECT

         ! setup transformation matrix to get an orthogonal basis, remove linear dependencies
         NULLIFY (integrals%utrans, integrals%uptrans)
         n = MAXVAL(basis%nbas)
         ALLOCATE (integrals%utrans(n, n, 0:lmat), integrals%uptrans(n, n, 0:lmat))
         integrals%utrans = 0._dp
         integrals%uptrans = 0._dp
         integrals%nne = integrals%n
         lwork = 10*n
         ALLOCATE (omat(n, n), vmat(n, n), w(n), work(lwork))
         DO l = 0, lmat
            n = integrals%n(l)
            IF (n > 0) THEN
               omat(1:n, 1:n) = integrals%ovlp(1:n, 1:n, l)
               CALL jacobi(omat(1:n, 1:n), w(1:n), vmat(1:n, 1:n))
               omat(1:n, 1:n) = vmat(1:n, 1:n)
               ii = 0
               DO i = 1, n
                  IF (w(i) > basis%eps_eig) THEN
                     ii = ii + 1
                     integrals%utrans(1:n, ii, l) = omat(1:n, i)/SQRT(w(i))
                  END IF
               END DO
               integrals%nne(l) = ii
               IF (ii > 0) THEN
                  omat(1:ii, 1:ii) = MATMUL(TRANSPOSE(integrals%utrans(1:n, 1:ii, l)), integrals%utrans(1:n, 1:ii, l))
                  DO i = 1, ii
                     integrals%uptrans(i, i, l) = 1._dp
                  ENDDO
                  CALL lapack_sgesv(ii, ii, omat(1:ii, 1:ii), ii, ipiv, integrals%uptrans(1:ii, 1:ii, l), ii, info)
                  CPASSERT(info == 0)
               END IF
            END IF
         END DO
         DEALLOCATE (omat, vmat, w, work)
      END IF

      CALL timestop(handle)

   END SUBROUTINE atom_int_setup

! **************************************************************************************************
!> \brief ...
!> \param integrals ...
!> \param basis ...
!> \param potential ...
! **************************************************************************************************
   SUBROUTINE atom_ppint_setup(integrals, basis, potential)
      TYPE(atom_integrals), INTENT(INOUT)                :: integrals
      TYPE(atom_basis_type), INTENT(INOUT)               :: basis
      TYPE(atom_potential_type), INTENT(IN)              :: potential

      CHARACTER(len=*), PARAMETER                        :: routineN = 'atom_ppint_setup'

      INTEGER                                            :: handle, i, ii, j, k, l, m, n
      REAL(KIND=dp)                                      :: al, alpha, rc
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: cpot, xmat
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: omat, spmat
      REAL(KIND=dp), DIMENSION(:), POINTER               :: rad

      CALL timeset(routineN, handle)

      IF (integrals%ppstat == 0) THEN
         n = MAXVAL(basis%nbas)
         integrals%n = basis%nbas

         NULLIFY (integrals%core, integrals%hnl)

         ALLOCATE (integrals%hnl(n, n, 0:lmat))
         integrals%hnl = 0._dp

         ALLOCATE (integrals%core(n, n, 0:lmat))
         integrals%core = 0._dp

         ALLOCATE (integrals%clsd(n, n, 0:lmat))
         integrals%clsd = 0._dp

         integrals%ppstat = 1

         SELECT CASE (basis%basis_type)
         CASE DEFAULT
            CPABORT("")
         CASE (GTO_BASIS)

            SELECT CASE (potential%ppot_type)
            CASE (no_pseudo, ecp_pseudo)
               DO l = 0, lmat
                  n = integrals%n(l)
                  CALL sg_nuclear(integrals%core(1:n, 1:n, l), l, basis%am(1:n, l), basis%am(1:n, l))
               END DO
            CASE (gth_pseudo)
               alpha = 1._dp/potential%gth_pot%rc/SQRT(2._dp)
               DO l = 0, lmat
                  n = integrals%n(l)
                  ALLOCATE (omat(n, n), spmat(n, 5))

                  omat = 0._dp
                  CALL sg_erf(omat(1:n, 1:n), l, alpha, basis%am(1:n, l), basis%am(1:n, l))
                  integrals%core(1:n, 1:n, l) = -potential%gth_pot%zion*omat(1:n, 1:n)
                  DO i = 1, potential%gth_pot%ncl
                     omat = 0._dp
                     CALL sg_gpot(omat(1:n, 1:n), i - 1, potential%gth_pot%rc, l, basis%am(1:n, l), basis%am(1:n, l))
                     integrals%core(1:n, 1:n, l) = integrals%core(1:n, 1:n, l) + &
                                                   potential%gth_pot%cl(i)*omat(1:n, 1:n)
                  END DO
                  IF (potential%gth_pot%lpotextended) THEN
                     DO k = 1, potential%gth_pot%nexp_lpot
                        DO i = 1, potential%gth_pot%nct_lpot(k)
                           omat = 0._dp
                           CALL sg_gpot(omat(1:n, 1:n), i - 1, potential%gth_pot%alpha_lpot(k), l, &
                                        basis%am(1:n, l), basis%am(1:n, l))
                           integrals%core(1:n, 1:n, l) = integrals%core(1:n, 1:n, l) + &
                                                         potential%gth_pot%cval_lpot(i, k)*omat(1:n, 1:n)
                        END DO
                     END DO
                  END IF
                  IF (potential%gth_pot%lsdpot) THEN
                     DO k = 1, potential%gth_pot%nexp_lsd
                        DO i = 1, potential%gth_pot%nct_lsd(k)
                           omat = 0._dp
                           CALL sg_gpot(omat(1:n, 1:n), i - 1, potential%gth_pot%alpha_lsd(k), l, &
                                        basis%am(1:n, l), basis%am(1:n, l))
                           integrals%clsd(1:n, 1:n, l) = integrals%clsd(1:n, 1:n, l) + &
                                                         potential%gth_pot%cval_lsd(i, k)*omat(1:n, 1:n)
                        END DO
                     END DO
                  END IF

                  spmat = 0._dp
                  m = potential%gth_pot%nl(l)
                  DO i = 1, m
                     CALL sg_proj_ol(spmat(1:n, i), l, basis%am(1:n, l), i - 1, potential%gth_pot%rcnl(l))
                  END DO
                  integrals%hnl(1:n, 1:n, l) = MATMUL(spmat(1:n, 1:m), &
                                                      MATMUL(potential%gth_pot%hnl(1:m, 1:m, l), TRANSPOSE(spmat(1:n, 1:m))))

                  DEALLOCATE (omat, spmat)
               END DO
            CASE (upf_pseudo)
               CALL upfint_setup(integrals, basis, potential)
            CASE (sgp_pseudo)
               CALL sgpint_setup(integrals, basis, potential)
            CASE DEFAULT
               CPABORT("")
            END SELECT

         CASE (CGTO_BASIS)

            SELECT CASE (potential%ppot_type)
            CASE (no_pseudo, ecp_pseudo)
               DO l = 0, lmat
                  n = integrals%n(l)
                  m = basis%nprim(l)
                  ALLOCATE (omat(m, m))

                  CALL sg_nuclear(omat(1:m, 1:m), l, basis%am(1:m, l), basis%am(1:m, l))
                  CALL contract2(integrals%core(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l))

                  DEALLOCATE (omat)
               END DO
            CASE (gth_pseudo)
               alpha = 1._dp/potential%gth_pot%rc/SQRT(2._dp)
               DO l = 0, lmat
                  n = integrals%n(l)
                  m = basis%nprim(l)
                  IF (n > 0 .AND. m > 0) THEN
                     ALLOCATE (omat(m, m), spmat(n, 5), xmat(m))

                     omat = 0._dp
                     CALL sg_erf(omat(1:m, 1:m), l, alpha, basis%am(1:m, l), basis%am(1:m, l))
                     omat(1:m, 1:m) = -potential%gth_pot%zion*omat(1:m, 1:m)
                     CALL contract2(integrals%core(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l))
                     DO i = 1, potential%gth_pot%ncl
                        omat = 0._dp
                        CALL sg_gpot(omat(1:m, 1:m), i - 1, potential%gth_pot%rc, l, basis%am(1:m, l), basis%am(1:m, l))
                        omat(1:m, 1:m) = potential%gth_pot%cl(i)*omat(1:m, 1:m)
                        CALL contract2add(integrals%core(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l))
                     END DO
                     IF (potential%gth_pot%lpotextended) THEN
                        DO k = 1, potential%gth_pot%nexp_lpot
                           DO i = 1, potential%gth_pot%nct_lpot(k)
                              omat = 0._dp
                              CALL sg_gpot(omat(1:m, 1:m), i - 1, potential%gth_pot%alpha_lpot(k), l, &
                                           basis%am(1:m, l), basis%am(1:m, l))
                              omat(1:m, 1:m) = potential%gth_pot%cval_lpot(i, k)*omat(1:m, 1:m)
                              CALL contract2add(integrals%core(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l))
                           END DO
                        END DO
                     END IF
                     IF (potential%gth_pot%lsdpot) THEN
                        DO k = 1, potential%gth_pot%nexp_lsd
                           DO i = 1, potential%gth_pot%nct_lsd(k)
                              omat = 0._dp
                              CALL sg_gpot(omat(1:m, 1:m), i - 1, potential%gth_pot%alpha_lsd(k), l, &
                                           basis%am(1:m, l), basis%am(1:m, l))
                              omat(1:m, 1:m) = potential%gth_pot%cval_lsd(i, k)*omat(1:m, 1:m)
                              CALL contract2add(integrals%clsd(1:n, 1:n, l), omat(1:m, 1:m), basis%cm(1:m, 1:n, l))
                           END DO
                        END DO
                     END IF

                     spmat = 0._dp
                     k = potential%gth_pot%nl(l)
                     DO i = 1, k
                        CALL sg_proj_ol(xmat(1:m), l, basis%am(1:m, l), i - 1, potential%gth_pot%rcnl(l))
                        spmat(1:n, i) = MATMUL(TRANSPOSE(basis%cm(1:m, 1:n, l)), xmat(1:m))
                     END DO
                     IF (k > 0) THEN
                        integrals%hnl(1:n, 1:n, l) = MATMUL(spmat(1:n, 1:k), &
                                                            MATMUL(potential%gth_pot%hnl(1:k, 1:k, l), &
                                                                   TRANSPOSE(spmat(1:n, 1:k))))
                     END IF
                     DEALLOCATE (omat, spmat, xmat)
                  END IF
               END DO
            CASE (upf_pseudo)
               CALL upfint_setup(integrals, basis, potential)
            CASE (sgp_pseudo)
               CALL sgpint_setup(integrals, basis, potential)
            CASE DEFAULT
               CPABORT("")
            END SELECT

         CASE (STO_BASIS)

            SELECT CASE (potential%ppot_type)
            CASE (no_pseudo, ecp_pseudo)
               DO l = 0, lmat
                  n = integrals%n(l)
                  CALL sto_nuclear(integrals%core(1:n, 1:n, l), basis%ns(1:n, l), basis%as(1:n, l), &
                                   basis%ns(1:n, l), basis%as(1:n, l))
               END DO
            CASE (gth_pseudo)
               rad => basis%grid%rad
               m = basis%grid%nr
               ALLOCATE (cpot(1:m))
               rc = potential%gth_pot%rc
               alpha = 1._dp/rc/SQRT(2._dp)

               ! local pseudopotential, we use erf = 1/r - erfc
               integrals%core = 0._dp
               DO i = 1, m
                  cpot(i) = potential%gth_pot%zion*erfc(alpha*rad(i))/rad(i)
               END DO
               DO i = 1, potential%gth_pot%ncl
                  ii = 2*(i - 1)
                  cpot(1:m) = cpot(1:m) + potential%gth_pot%cl(i)*(rad/rc)**ii*EXP(-0.5_dp*(rad/rc)**2)
               END DO
               IF (potential%gth_pot%lpotextended) THEN
                  DO k = 1, potential%gth_pot%nexp_lpot
                     al = potential%gth_pot%alpha_lpot(k)
                     DO i = 1, potential%gth_pot%nct_lpot(k)
                        ii = 2*(i - 1)
                        cpot(1:m) = cpot(1:m) + potential%gth_pot%cval_lpot(i, k)*(rad/al)**ii*EXP(-0.5_dp*(rad/al)**2)
                     END DO
                  END DO
               END IF
               CALL numpot_matrix(integrals%core, cpot, basis, 0)
               DO l = 0, lmat
                  n = integrals%n(l)
                  ALLOCATE (omat(n, n))
                  omat = 0._dp
                  CALL sto_nuclear(omat(1:n, 1:n), basis%ns(1:n, l), basis%as(1:n, l), &
                                   basis%ns(1:n, l), basis%as(1:n, l))
                  integrals%core(1:n, 1:n, l) = integrals%core(1:n, 1:n, l) - potential%gth_pot%zion*omat(1:n, 1:n)
                  DEALLOCATE (omat)
               END DO

               IF (potential%gth_pot%lsdpot) THEN
                  cpot = 0._dp
                  DO k = 1, potential%gth_pot%nexp_lsd
                     al = potential%gth_pot%alpha_lsd(k)
                     DO i = 1, potential%gth_pot%nct_lsd(k)
                        ii = 2*(i - 1)
                        cpot(:) = cpot + potential%gth_pot%cval_lsd(i, k)*(rad/al)**ii*EXP(-0.5_dp*(rad/al)**2)
                     END DO
                  END DO
                  CALL numpot_matrix(integrals%clsd, cpot, basis, 0)
               END IF

               DO l = 0, lmat
                  n = integrals%n(l)
                  ! non local pseudopotential
                  ALLOCATE (spmat(n, 5))
                  spmat = 0._dp
                  k = potential%gth_pot%nl(l)
                  DO i = 1, k
                     rc = potential%gth_pot%rcnl(l)
                     cpot(:) = sqrt2/SQRT(gamma1(l + 2*i - 1))*rad**(l + 2*i - 2)*EXP(-0.5_dp*(rad/rc)**2)/rc**(l + 2*i - 0.5_dp)
                     DO j = 1, basis%nbas(l)
                        spmat(j, i) = integrate_grid(cpot, basis%bf(:, j, l), basis%grid)
                     END DO
                  END DO
                  integrals%hnl(1:n, 1:n, l) = MATMUL(spmat(1:n, 1:k), &
                                                      MATMUL(potential%gth_pot%hnl(1:k, 1:k, l), &
                                                             TRANSPOSE(spmat(1:n, 1:k))))
                  DEALLOCATE (spmat)
               END DO
               DEALLOCATE (cpot)

            CASE (upf_pseudo)
               CALL upfint_setup(integrals, basis, potential)
            CASE (sgp_pseudo)
               CALL sgpint_setup(integrals, basis, potential)
            CASE DEFAULT
               CPABORT("")
            END SELECT

         CASE (NUM_BASIS)
            CPABORT("")
         END SELECT

         ! add ecp_pseudo using numerical representation of basis
         IF (potential%ppot_type == ecp_pseudo) THEN
            ! scale 1/r potential
            integrals%core = -potential%ecp_pot%zion*integrals%core
            ! local potential
            m = basis%grid%nr
            rad => basis%grid%rad
            ALLOCATE (cpot(1:m))
            cpot = 0._dp
            DO k = 1, potential%ecp_pot%nloc
               n = potential%ecp_pot%nrloc(k)
               alpha = potential%ecp_pot%bloc(k)
               cpot(:) = cpot + potential%ecp_pot%aloc(k)*rad**(n - 2)*EXP(-alpha*rad**2)
            END DO
            CALL numpot_matrix(integrals%core, cpot, basis, 0)
            ! non local pseudopotential
            DO l = 0, MIN(potential%ecp_pot%lmax, lmat)
               cpot = 0._dp
               DO k = 1, potential%ecp_pot%npot(l)
                  n = potential%ecp_pot%nrpot(k, l)
                  alpha = potential%ecp_pot%bpot(k, l)
                  cpot(:) = cpot + potential%ecp_pot%apot(k, l)*rad**(n - 2)*EXP(-alpha*rad**2)
               END DO
               DO i = 1, basis%nbas(l)
                  DO j = i, basis%nbas(l)
                     integrals%hnl(i, j, l) = integrate_grid(cpot, basis%bf(:, i, l), basis%bf(:, j, l), basis%grid)
                     integrals%hnl(j, i, l) = integrals%hnl(i, j, l)
                  END DO
               END DO
            END DO
            DEALLOCATE (cpot)
         END IF

      END IF

      CALL timestop(handle)

   END SUBROUTINE atom_ppint_setup

! **************************************************************************************************
!> \brief ...
!> \param integrals ...
!> \param basis ...
!> \param potential ...
! **************************************************************************************************
   SUBROUTINE upfint_setup(integrals, basis, potential)
      TYPE(atom_integrals), INTENT(INOUT)                :: integrals
      TYPE(atom_basis_type), INTENT(INOUT)               :: basis
      TYPE(atom_potential_type), INTENT(IN)              :: potential

      CHARACTER(len=4)                                   :: ptype
      INTEGER                                            :: i, j, k1, k2, la, lb, m, n
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: spot
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: spmat
      TYPE(atom_basis_type)                              :: gbasis

      ! get basis representation on UPF grid
      CALL atom_basis_gridrep(basis, gbasis, potential%upf_pot%r, potential%upf_pot%rab)

      ! local pseudopotential
      integrals%core = 0._dp
      CALL numpot_matrix(integrals%core, potential%upf_pot%vlocal, gbasis, 0)

      ptype = ADJUSTL(TRIM(potential%upf_pot%pseudo_type))
      IF (ptype(1:2) == "NC" .OR. ptype(1:2) == "US") THEN
         ! non local pseudopotential
         n = MAXVAL(integrals%n(:))
         m = potential%upf_pot%number_of_proj
         ALLOCATE (spmat(n, m))
         spmat = 0.0_dp
         DO i = 1, m
            la = potential%upf_pot%lbeta(i)
            DO j = 1, gbasis%nbas(la)
               spmat(j, i) = integrate_grid(potential%upf_pot%beta(:, i), gbasis%bf(:, j, la), gbasis%grid)
            END DO
         END DO
         DO i = 1, m
            la = potential%upf_pot%lbeta(i)
            DO j = 1, m
               lb = potential%upf_pot%lbeta(j)
               IF (la == lb) THEN
                  DO k1 = 1, gbasis%nbas(la)
                     DO k2 = 1, gbasis%nbas(la)
                        integrals%hnl(k1, k2, la) = integrals%hnl(k1, k2, la) + &
                                                    spmat(k1, i)*potential%upf_pot%dion(i, j)*spmat(k2, j)
                     END DO
                  END DO
               END IF
            END DO
         END DO
         DEALLOCATE (spmat)
      ELSE IF (ptype(1:2) == "SL") THEN
         ! semi local pseudopotential
         DO la = 0, potential%upf_pot%l_max
            IF (la == potential%upf_pot%l_local) CYCLE
            m = SIZE(potential%upf_pot%vsemi(:, la + 1))
            ALLOCATE (spot(m))
            spot(:) = potential%upf_pot%vsemi(:, la + 1) - potential%upf_pot%vlocal(:)
            n = basis%nbas(la)
            DO i = 1, n
               DO j = i, n
                  integrals%core(i, j, la) = integrals%core(i, j, la) + &
                                             integrate_grid(spot(:), &
                                                            gbasis%bf(:, i, la), gbasis%bf(:, j, la), gbasis%grid)
                  integrals%core(j, i, la) = integrals%core(i, j, la)
               END DO
            END DO
            DEALLOCATE (spot)
         END DO
      ELSE
         CPABORT("Pseudopotential type: ["//ADJUSTL(TRIM(ptype))//"] not known")
      END IF

      ! release basis representation on UPF grid
      CALL release_atom_basis(gbasis)

   END SUBROUTINE upfint_setup

! **************************************************************************************************
!> \brief ...
!> \param integrals ...
!> \param basis ...
!> \param potential ...
! **************************************************************************************************
   SUBROUTINE sgpint_setup(integrals, basis, potential)
      TYPE(atom_integrals), INTENT(INOUT)                :: integrals
      TYPE(atom_basis_type), INTENT(INOUT)               :: basis
      TYPE(atom_potential_type), INTENT(IN)              :: potential

      INTEGER                                            :: i, ia, j, l, m, n, na
      REAL(KIND=dp)                                      :: a, c, rc, zval
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: cpot, pgauss
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: qmat
      REAL(KIND=dp), DIMENSION(:), POINTER               :: rad

      rad => basis%grid%rad
      m = basis%grid%nr

      ! local pseudopotential
      integrals%core = 0._dp
      ALLOCATE (cpot(m))
      cpot = 0.0_dp
      zval = potential%sgp_pot%zion
      DO i = 1, m
         rc = rad(i)/potential%sgp_pot%ac_local/SQRT(2.0_dp)
         cpot(i) = cpot(i) - zval/rad(i)*erf(rc)
      END DO
      DO i = 1, potential%sgp_pot%n_local
         cpot(:) = cpot(:) + potential%sgp_pot%c_local(i)*EXP(-potential%sgp_pot%a_local(i)*rad(:)**2)
      END DO
      CALL numpot_matrix(integrals%core, cpot, basis, 0)
      DEALLOCATE (cpot)

      ! nonlocal pseudopotential
      integrals%hnl = 0.0_dp
      IF (potential%sgp_pot%has_nonlocal) THEN
         ALLOCATE (pgauss(1:m))
         n = potential%sgp_pot%n_nonlocal
         !
         DO l = 0, potential%sgp_pot%lmax
            CPASSERT(l <= UBOUND(basis%nbas, 1))
            IF (.NOT. potential%sgp_pot%is_nonlocal(l)) CYCLE
            ! overlap (a|p)
            na = basis%nbas(l)
            ALLOCATE (qmat(na, n))
            DO i = 1, n
               pgauss(:) = 0.0_dp
               DO j = 1, n
                  a = potential%sgp_pot%a_nonlocal(j)
                  c = potential%sgp_pot%c_nonlocal(j, i, l)
                  pgauss(:) = pgauss(:) + c*EXP(-a*rad(:)**2)*rad(:)**l
               END DO
               DO ia = 1, na
                  qmat(ia, i) = SUM(basis%bf(:, ia, l)*pgauss(:)*basis%grid%wr(:))
               END DO
            END DO
            DO i = 1, na
               DO j = i, na
                  DO ia = 1, n
                     integrals%hnl(i, j, l) = integrals%hnl(i, j, l) &
                                              + qmat(i, ia)*qmat(j, ia)*potential%sgp_pot%h_nonlocal(ia, l)
                  END DO
                  integrals%hnl(j, i, l) = integrals%hnl(i, j, l)
               END DO
            END DO
            DEALLOCATE (qmat)
         END DO
         DEALLOCATE (pgauss)
      END IF

   END SUBROUTINE sgpint_setup

! **************************************************************************************************
!> \brief ...
!> \param integrals ...
!> \param basis ...
!> \param reltyp ...
!> \param zcore ...
!> \param alpha ...
! **************************************************************************************************
   SUBROUTINE atom_relint_setup(integrals, basis, reltyp, zcore, alpha)
      TYPE(atom_integrals), INTENT(INOUT)                :: integrals
      TYPE(atom_basis_type), INTENT(INOUT)               :: basis
      INTEGER, INTENT(IN)                                :: reltyp
      REAL(dp), INTENT(IN)                               :: zcore
      REAL(dp), INTENT(IN), OPTIONAL                     :: alpha

      CHARACTER(len=*), PARAMETER                        :: routineN = 'atom_relint_setup'

      INTEGER                                            :: dkhorder, handle, i, k1, k2, l, m, n, nl
      REAL(dp)                                           :: ascal
      REAL(dp), ALLOCATABLE, DIMENSION(:)                :: cpot
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: modpot
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: ener, sps
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: hmat, pvp, sp, tp, vp, wfn

      CALL timeset(routineN, handle)

      SELECT CASE (reltyp)
      CASE DEFAULT
         CPABORT("")
      CASE (do_nonrel_atom, do_zoramp_atom, do_sczoramp_atom)
         dkhorder = -1
      CASE (do_dkh0_atom)
         dkhorder = 0
      CASE (do_dkh1_atom)
         dkhorder = 1
      CASE (do_dkh2_atom)
         dkhorder = 2
      CASE (do_dkh3_atom)
         dkhorder = 3
      END SELECT

      SELECT CASE (reltyp)
      CASE DEFAULT
         CPABORT("")
      CASE (do_nonrel_atom)
         ! nothing to do
         NULLIFY (integrals%tzora, integrals%hdkh)
      CASE (do_zoramp_atom, do_sczoramp_atom)

         NULLIFY (integrals%hdkh)

         IF (integrals%zorastat == 0) THEN
            n = MAXVAL(basis%nbas)
            ALLOCATE (integrals%tzora(n, n, 0:lmat))
            integrals%tzora = 0._dp
            m = basis%grid%nr
            ALLOCATE (modpot(1:m), cpot(1:m))
            CALL calculate_model_potential(modpot, basis%grid, zcore)
            ! Zora potential
            cpot(1:m) = modpot(1:m)/(4._dp*c_light_au*c_light_au - 2._dp*modpot(1:m))
            cpot(1:m) = cpot(1:m)/basis%grid%rad2(1:m)
            CALL numpot_matrix(integrals%tzora, cpot, basis, 0)
            DO l = 0, lmat
               nl = basis%nbas(l)
               integrals%tzora(1:nl, 1:nl, l) = REAL(l*(l + 1), dp)*integrals%tzora(1:nl, 1:nl, l)
            END DO
            cpot(1:m) = cpot(1:m)*basis%grid%rad2(1:m)
            CALL numpot_matrix(integrals%tzora, cpot, basis, 2)
            !
            ! scaled ZORA
            IF (reltyp == do_sczoramp_atom) THEN
               ALLOCATE (hmat(n, n, 0:lmat), wfn(n, n, 0:lmat), ener(n, 0:lmat), pvp(n, n, 0:lmat), sps(n, n))
               hmat(:, :, :) = integrals%kin + integrals%tzora
               ! model potential
               CALL numpot_matrix(hmat, modpot, basis, 0)
               ! eigenvalues and eigenvectors
               CALL atom_solve(hmat, integrals%utrans, wfn, ener, basis%nbas, integrals%nne, lmat)
               ! relativistic kinetic energy
               cpot(1:m) = c_light_au*c_light_au/(2._dp*c_light_au*c_light_au - modpot(1:m))**2
               cpot(1:m) = cpot(1:m)/basis%grid%rad2(1:m)
               pvp = 0.0_dp
               CALL numpot_matrix(pvp, cpot, basis, 0)
               DO l = 0, lmat
                  nl = basis%nbas(l)
                  pvp(1:nl, 1:nl, l) = REAL(l*(l + 1), dp)*pvp(1:nl, 1:nl, l)
               END DO
               cpot(1:m) = cpot(1:m)*basis%grid%rad2(1:m)
               CALL numpot_matrix(pvp, cpot, basis, 2)
               ! calculate psi*pvp*psi and the scaled orbital energies
               ! actually, we directly calculate the energy difference
               DO l = 0, lmat
                  nl = basis%nbas(l)
                  DO i = 1, integrals%nne(l)
                     IF (ener(i, l) < 0._dp) THEN
                        ascal = SUM(wfn(1:nl, i, l)*MATMUL(pvp(1:nl, 1:nl, l), wfn(1:nl, i, l)))
                        ener(i, l) = ener(i, l)*ascal/(1.0_dp + ascal)
                     ELSE
                        ener(i, l) = 0.0_dp
                     END IF
                  END DO
               END DO
               ! correction term is calculated as a projector
               hmat = 0.0_dp
               DO l = 0, lmat
                  nl = basis%nbas(l)
                  DO i = 1, integrals%nne(l)
                     DO k1 = 1, nl
                        DO k2 = 1, nl
                           hmat(k1, k2, l) = hmat(k1, k2, l) + ener(i, l)*wfn(k1, i, l)*wfn(k2, i, l)
                        END DO
                     END DO
                  END DO
                  ! transform with overlap matrix
                  sps(1:nl, 1:nl) = MATMUL(integrals%ovlp(1:nl, 1:nl, l), &
                                           MATMUL(hmat(1:nl, 1:nl, l), integrals%ovlp(1:nl, 1:nl, l)))
                  ! add scaling correction to tzora
                  integrals%tzora(1:nl, 1:nl, l) = integrals%tzora(1:nl, 1:nl, l) - sps(1:nl, 1:nl)
               END DO

               DEALLOCATE (hmat, wfn, ener, pvp, sps)
            END IF
            !
            DEALLOCATE (modpot, cpot)

            integrals%zorastat = 1

         END IF

      CASE (do_dkh0_atom, do_dkh1_atom, do_dkh2_atom, do_dkh3_atom)

         NULLIFY (integrals%tzora)

         IF (integrals%dkhstat == 0) THEN
            n = MAXVAL(basis%nbas)
            ALLOCATE (integrals%hdkh(n, n, 0:lmat))
            integrals%hdkh = 0._dp

            m = MAXVAL(basis%nprim)
            ALLOCATE (tp(m, m, 0:lmat), sp(m, m, 0:lmat), vp(m, m, 0:lmat), pvp(m, m, 0:lmat))
            tp = 0._dp
            sp = 0._dp
            vp = 0._dp
            pvp = 0._dp

            SELECT CASE (basis%basis_type)
            CASE DEFAULT
               CPABORT("")
            CASE (GTO_BASIS, CGTO_BASIS)

               DO l = 0, lmat
                  m = basis%nprim(l)
                  IF (m > 0) THEN
                     CALL sg_kinetic(tp(1:m, 1:m, l), l, basis%am(1:m, l), basis%am(1:m, l))
                     CALL sg_overlap(sp(1:m, 1:m, l), l, basis%am(1:m, l), basis%am(1:m, l))
                     IF (PRESENT(alpha)) THEN
                        CALL sg_erfc(vp(1:m, 1:m, l), l, alpha, basis%am(1:m, l), basis%am(1:m, l))
                     ELSE
                        CALL sg_nuclear(vp(1:m, 1:m, l), l, basis%am(1:m, l), basis%am(1:m, l))
                     END IF
                     CALL sg_kinnuc(pvp(1:m, 1:m, l), l, basis%am(1:m, l), basis%am(1:m, l))
                     vp(1:m, 1:m, l) = -zcore*vp(1:m, 1:m, l)
                     pvp(1:m, 1:m, l) = -zcore*pvp(1:m, 1:m, l)
                  END IF
               END DO

            CASE (STO_BASIS)
               CPABORT("")
            CASE (NUM_BASIS)
               CPABORT("")
            END SELECT

            CALL dkh_integrals(integrals, basis, dkhorder, sp, tp, vp, pvp)

            integrals%dkhstat = 1

            DEALLOCATE (tp, sp, vp, pvp)

         ELSE
            CPASSERT(ASSOCIATED(integrals%hdkh))
         END IF

      END SELECT

      CALL timestop(handle)

   END SUBROUTINE atom_relint_setup

! **************************************************************************************************
!> \brief ...
!> \param integrals ...
!> \param basis ...
!> \param order ...
!> \param sp ...
!> \param tp ...
!> \param vp ...
!> \param pvp ...
! **************************************************************************************************
   SUBROUTINE dkh_integrals(integrals, basis, order, sp, tp, vp, pvp)
      TYPE(atom_integrals), INTENT(INOUT)                :: integrals
      TYPE(atom_basis_type), INTENT(INOUT)               :: basis
      INTEGER, INTENT(IN)                                :: order
      REAL(dp), DIMENSION(:, :, 0:)                      :: sp, tp, vp, pvp

      INTEGER                                            :: l, m, n
      REAL(dp), DIMENSION(:, :, :), POINTER              :: hdkh

      CPASSERT(order >= 0)

      hdkh => integrals%hdkh

      DO l = 0, lmat
         n = integrals%n(l)
         m = basis%nprim(l)
         IF (n > 0) THEN
            CALL dkh_atom_transformation(sp(1:m, 1:m, l), vp(1:m, 1:m, l), tp(1:m, 1:m, l), pvp(1:m, 1:m, l), m, order)
            SELECT CASE (basis%basis_type)
            CASE DEFAULT
               CPABORT("")
            CASE (GTO_BASIS)
               CPASSERT(n == m)
               integrals%hdkh(1:n, 1:n, l) = tp(1:n, 1:n, l) + vp(1:n, 1:n, l)
            CASE (CGTO_BASIS)
               CALL contract2(integrals%hdkh(1:n, 1:n, l), tp(1:m, 1:m, l), basis%cm(1:m, 1:n, l))
               CALL contract2add(integrals%hdkh(1:n, 1:n, l), vp(1:m, 1:m, l), basis%cm(1:m, 1:n, l))
            CASE (STO_BASIS)
               CPABORT("")
            CASE (NUM_BASIS)
               CPABORT("")
            END SELECT
         ELSE
            integrals%hdkh(1:n, 1:n, l) = 0._dp
         END IF
      END DO

   END SUBROUTINE dkh_integrals

! **************************************************************************************************
!> \brief Calculate overlap matrix between two atomic basis sets.
!> \param ovlap    overlap matrix <chi_{a,l} | chi_{b,l}>
!> \param basis_a  first basis set
!> \param basis_b  second basis set
! **************************************************************************************************
   SUBROUTINE atom_basis_projection_overlap(ovlap, basis_a, basis_b)
      REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(OUT)    :: ovlap
      TYPE(atom_basis_type), INTENT(IN)                  :: basis_a, basis_b

      INTEGER                                            :: i, j, l, ma, mb, na, nb
      LOGICAL                                            :: ebas
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: omat

      ebas = (basis_a%basis_type == basis_b%basis_type)

      ovlap = 0.0_dp

      IF (ebas) THEN
         SELECT CASE (basis_a%basis_type)
         CASE DEFAULT
            CPABORT("")
         CASE (GTO_BASIS)
            DO l = 0, lmat
               na = basis_a%nbas(l)
               nb = basis_b%nbas(l)
               CALL sg_overlap(ovlap(1:na, 1:nb, l), l, basis_a%am(1:na, l), basis_b%am(1:nb, l))
            END DO
         CASE (CGTO_BASIS)
            DO l = 0, lmat
               na = basis_a%nbas(l)
               nb = basis_b%nbas(l)
               ma = basis_a%nprim(l)
               mb = basis_b%nprim(l)
               ALLOCATE (omat(ma, mb))
               CALL sg_overlap(omat(1:ma, 1:mb), l, basis_a%am(1:ma, l), basis_b%am(1:mb, l))
               ovlap(1:na, 1:nb, l) = MATMUL(TRANSPOSE(basis_a%cm(1:ma, 1:na, l)), &
                                             MATMUL(omat(1:ma, 1:mb), basis_b%cm(1:mb, 1:nb, l)))
               DEALLOCATE (omat)
            END DO
         CASE (STO_BASIS)
            DO l = 0, lmat
               na = basis_a%nbas(l)
               nb = basis_b%nbas(l)
               CALL sto_overlap(ovlap(1:na, 1:nb, l), basis_a%ns(1:na, l), basis_b%as(1:nb, l), &
                                basis_a%ns(1:na, l), basis_b%as(1:nb, l))
            END DO
         CASE (NUM_BASIS)
            CPASSERT(atom_compare_grids(basis_a%grid, basis_b%grid))
            DO l = 0, lmat
               na = basis_a%nbas(l)
               nb = basis_b%nbas(l)
               DO j = 1, nb
                  DO i = 1, na
                     ovlap(i, j, l) = integrate_grid(basis_a%bf(:, i, l), basis_b%bf(:, j, l), basis_a%grid)
                  END DO
               END DO
            END DO
         END SELECT
      ELSE
         CPASSERT(atom_compare_grids(basis_a%grid, basis_b%grid))
         DO l = 0, lmat
            na = basis_a%nbas(l)
            nb = basis_b%nbas(l)
            DO j = 1, nb
               DO i = 1, na
                  ovlap(i, j, l) = integrate_grid(basis_a%bf(:, i, l), basis_b%bf(:, j, l), basis_a%grid)
               END DO
            END DO
         END DO
      END IF

   END SUBROUTINE atom_basis_projection_overlap

! **************************************************************************************************
!> \brief Release memory allocated for atomic integrals (valence electrons).
!> \param integrals  atomic integrals
! **************************************************************************************************
   SUBROUTINE atom_int_release(integrals)
      TYPE(atom_integrals), INTENT(INOUT)                :: integrals

      INTEGER                                            :: ll

      IF (ASSOCIATED(integrals%ovlp)) THEN
         DEALLOCATE (integrals%ovlp)
      END IF
      IF (ASSOCIATED(integrals%kin)) THEN
         DEALLOCATE (integrals%kin)
      END IF
      IF (ASSOCIATED(integrals%conf)) THEN
         DEALLOCATE (integrals%conf)
      END IF
      DO ll = 1, SIZE(integrals%ceri)
         IF (ASSOCIATED(integrals%ceri(ll)%int)) THEN
            DEALLOCATE (integrals%ceri(ll)%int)
         END IF
         IF (ASSOCIATED(integrals%eeri(ll)%int)) THEN
            DEALLOCATE (integrals%eeri(ll)%int)
         END IF
      END DO
      IF (ASSOCIATED(integrals%utrans)) THEN
         DEALLOCATE (integrals%utrans)
      END IF
      IF (ASSOCIATED(integrals%uptrans)) THEN
         DEALLOCATE (integrals%uptrans)
      END IF

      integrals%status = 0

   END SUBROUTINE atom_int_release

! **************************************************************************************************
!> \brief Release memory allocated for atomic integrals (core electrons).
!> \param integrals  atomic integrals
! **************************************************************************************************
   SUBROUTINE atom_ppint_release(integrals)
      TYPE(atom_integrals), INTENT(INOUT)                :: integrals

      IF (ASSOCIATED(integrals%hnl)) THEN
         DEALLOCATE (integrals%hnl)
      END IF
      IF (ASSOCIATED(integrals%core)) THEN
         DEALLOCATE (integrals%core)
      END IF
      IF (ASSOCIATED(integrals%clsd)) THEN
         DEALLOCATE (integrals%clsd)
      END IF

      integrals%ppstat = 0

   END SUBROUTINE atom_ppint_release

! **************************************************************************************************
!> \brief  Release memory allocated for atomic integrals (relativistic effects).
!> \param integrals atomic integrals
! **************************************************************************************************
   SUBROUTINE atom_relint_release(integrals)
      TYPE(atom_integrals), INTENT(INOUT)                :: integrals

      IF (ASSOCIATED(integrals%tzora)) THEN
         DEALLOCATE (integrals%tzora)
      END IF
      IF (ASSOCIATED(integrals%hdkh)) THEN
         DEALLOCATE (integrals%hdkh)
      END IF

      integrals%zorastat = 0
      integrals%dkhstat = 0

   END SUBROUTINE atom_relint_release

! **************************************************************************************************
!> \brief Calculate model potential. It is useful to guess atomic orbitals.
!> \param modpot  model potential
!> \param grid    atomic radial grid
!> \param zcore   nuclear charge
! **************************************************************************************************
   SUBROUTINE calculate_model_potential(modpot, grid, zcore)
      REAL(dp), DIMENSION(:), INTENT(INOUT)              :: modpot
      TYPE(grid_atom_type), INTENT(IN)                   :: grid
      REAL(dp), INTENT(IN)                               :: zcore

      INTEGER                                            :: i, ii, l, ll, n, z
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: pot, rho
      TYPE(atom_state)                                   :: state

      n = SIZE(modpot)
      ALLOCATE (rho(n), pot(n))

      ! fill default occupation
      state%core = 0._dp
      state%occ = 0._dp
      state%multiplicity = -1
      z = NINT(zcore)
      DO l = 0, MIN(lmat, UBOUND(ptable(z)%e_conv, 1))
         IF (ptable(z)%e_conv(l) /= 0) THEN
            state%maxl_occ = l
            ll = 2*(2*l + 1)
            DO i = 1, SIZE(state%occ, 2)
               ii = ptable(z)%e_conv(l) - (i - 1)*ll
               IF (ii <= ll) THEN
                  state%occ(l, i) = ii
                  EXIT
               ELSE
                  state%occ(l, i) = ll
               END IF
            END DO
         END IF
      END DO

      modpot = -zcore/grid%rad(:)

      ! Coulomb potential
      CALL slater_density(rho, pot, NINT(zcore), state, grid)
      CALL coulomb_potential_numeric(pot, rho, grid)
      modpot = modpot + pot

      ! XC potential
      CALL wigner_slater_functional(rho, pot)
      modpot = modpot + pot

      DEALLOCATE (rho, pot)

   END SUBROUTINE calculate_model_potential

END MODULE atom_operators
