atom_sgp.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                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
MODULE atom_sgp
   USE atom_set_basis,                  ONLY: set_kind_basis_atomic
   USE atom_types,                      ONLY: &
        atom_basis_gridrep, atom_basis_type, atom_ecppot_type, atom_p_type, atom_type, &
        create_opmat, init_atom_basis_default_pp, opmat_type, release_atom_basis, release_opmat
   USE atom_utils,                      ONLY: integrate_grid,&
                                              numpot_matrix
   USE basis_set_types,                 ONLY: gto_basis_set_type
   USE input_constants,                 ONLY: ecp_pseudo,&
                                              gaussian,&
                                              gth_pseudo,&
                                              no_pseudo,&
                                              upf_pseudo
   USE input_section_types,             ONLY: section_vals_get,&
                                              section_vals_type
   USE kahan_sum,                       ONLY: accurate_dot_product
   USE kinds,                           ONLY: dp
   USE mathconstants,                   ONLY: dfac,&
                                              fourpi,&
                                              rootpi,&
                                              sqrt2
   USE mathlib,                         ONLY: diamat_all,&
                                              get_pseudo_inverse_diag
   USE powell,                          ONLY: opt_state_type,&
                                              powell_optimize
   USE upf_types,                       ONLY: atom_upfpot_type
#include "../base/base_uses.f90"
 
   IMPLICIT NONE
 
   TYPE atom_sgp_potential_type
      LOGICAL                                  :: has_nonlocal = .false.
      INTEGER                                  :: n_nonlocal = 0
      INTEGER                                  :: lmax = 0
      LOGICAL, DIMENSION(0:5)                  :: is_nonlocal = .false.
      REAL(kind=dp), DIMENSION(:), POINTER     :: a_nonlocal => null()
      REAL(kind=dp), DIMENSION(:, :), POINTER  :: h_nonlocal => null()
      REAL(kind=dp), DIMENSION(:, :, :), POINTER  :: c_nonlocal => null()
      LOGICAL                                  :: has_local = .false.
      INTEGER                                  :: n_local = 0
      REAL(kind=dp)                            :: zval = 0.0_dp
      REAL(kind=dp)                            :: ac_local = 0.0_dp
      REAL(kind=dp), DIMENSION(:), POINTER     :: a_local => null()
      REAL(kind=dp), DIMENSION(:), POINTER     :: c_local => null()
      LOGICAL                                  :: has_nlcc = .false.
      INTEGER                                  :: n_nlcc = 0
      REAL(kind=dp), DIMENSION(:), POINTER     :: a_nlcc => null()
      REAL(kind=dp), DIMENSION(:), POINTER     :: c_nlcc => null()
   END TYPE atom_sgp_potential_type
 
   PRIVATE
   PUBLIC  :: atom_sgp_potential_type, atom_sgp_release
   PUBLIC  :: atom_sgp_construction, sgp_construction
   PUBLIC  :: fit_core_density_gaussians
 
   CHARACTER(len=*), PARAMETER, PRIVATE :: modulen = 'atom_sgp'
 
! **************************************************************************************************
 
CONTAINS
 
! **************************************************************************************************
!> \brief ...
!> \param sgp_pot ...
!> \param ecp_pot ...
!> \param upf_pot ...
!> \param orb_basis ...
!> \param error ...
! **************************************************************************************************
   SUBROUTINE sgp_construction(sgp_pot, ecp_pot, upf_pot, orb_basis, error)
 
      TYPE(atom_sgp_potential_type)                      :: sgp_pot
      TYPE(atom_ecppot_type), OPTIONAL                   :: ecp_pot
      TYPE(atom_upfpot_type), OPTIONAL                   :: upf_pot
      TYPE(gto_basis_set_type), OPTIONAL, POINTER        :: orb_basis
      REAL(kind=dp), DIMENSION(6)                        :: error
 
      INTEGER                                            :: i, n
      INTEGER, DIMENSION(3)                              :: mloc
      LOGICAL                                            :: is_ecp, is_upf
      REAL(kind=dp)                                      :: errcc, rcut
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: cgauss, cutpots, cutpotu
      TYPE(atom_basis_type)                              :: basis
      TYPE(opmat_type), POINTER                          :: core, hnl, score, shnl
 
      ! define basis
      IF (PRESENT(orb_basis)) THEN
         CALL set_kind_basis_atomic(basis, orb_basis, has_pp=.true., cp2k_norm=.true.)
      ELSE
         CALL init_atom_basis_default_pp(basis)
      END IF
 
      is_ecp = .false.
      IF (PRESENT(ecp_pot)) is_ecp = .true.
      is_upf = .false.
      IF (PRESENT(upf_pot)) is_upf = .true.
      cpassert(.NOT. (is_ecp .AND. is_upf))
 
      ! upf has often very small grids, use a smooth cutoff function
      IF (is_upf) THEN
         n = SIZE(upf_pot%r)
         ALLOCATE (cutpotu(n))
         rcut = maxval(upf_pot%r)
         CALL cutpot(cutpotu, upf_pot%r, rcut, 2.5_dp)
         n = basis%grid%nr
         ALLOCATE (cutpots(n))
         CALL cutpot(cutpots, basis%grid%rad, rcut, 2.5_dp)
      ELSE
         n = basis%grid%nr
         ALLOCATE (cutpots(n))
         cutpots = 1.0_dp
      END IF
 
      ! generate the transformed potentials
      IF (is_ecp) THEN
         CALL ecp_sgp_constr(ecp_pot, sgp_pot, basis)
      ELSE IF (is_upf) THEN
         CALL upf_sgp_constr(upf_pot, sgp_pot, basis)
      ELSE
         cpabort("Either ecp_pot or upf_pot is needed for sgp_construction")
      END IF
 
      NULLIFY (core, hnl)
      CALL create_opmat(core, basis%nbas)
      CALL create_opmat(hnl, basis%nbas, 5)
      NULLIFY (score, shnl)
      CALL create_opmat(score, basis%nbas)
      CALL create_opmat(shnl, basis%nbas, 5)
      !
      IF (is_ecp) THEN
         CALL ecpints(hnl%op, basis, ecp_pot)
      ELSE IF (is_upf) THEN
         CALL upfints(core%op, hnl%op, basis, upf_pot, cutpotu, sgp_pot%ac_local)
      ELSE
         cpabort("Either ecp_pot or upf_pot is needed for sgp_construction")
      END IF
      !
      CALL sgpints(score%op, shnl%op, basis, sgp_pot, cutpots)
      !
      error = 0.0_dp
      IF (sgp_pot%has_local) THEN
         n = min(3, ubound(core%op, 3))
         error(1) = maxval(abs(core%op(:, :, 0:n) - score%op(:, :, 0:n)))
         mloc = maxloc(abs(core%op(:, :, 0:n) - score%op(:, :, 0:n)))
         ! MAXLOC returns positions starting at one, whereas the angular index starts at zero.
         mloc(3) = mloc(3) - 1
         error(4) = error(1)/max(abs(core%op(mloc(1), mloc(2), mloc(3))), 1.e-6_dp)
      END IF
      IF (sgp_pot%has_nonlocal) THEN
         n = min(3, ubound(hnl%op, 3))
         error(2) = maxval(abs(hnl%op(:, :, 0:n) - shnl%op(:, :, 0:n)))
         mloc = maxloc(abs(hnl%op(:, :, 0:n) - shnl%op(:, :, 0:n)))
         mloc(3) = mloc(3) - 1
         error(5) = error(2)/max(abs(hnl%op(mloc(1), mloc(2), mloc(3))), 1.e-6_dp)
      END IF
      IF (sgp_pot%has_nlcc) THEN
         IF (is_upf) THEN
            n = SIZE(upf_pot%r)
            ALLOCATE (cgauss(n))
            cgauss = 0.0_dp
            DO i = 1, sgp_pot%n_nlcc
               cgauss(:) = cgauss(:) + sgp_pot%c_nlcc(i)*exp(-sgp_pot%a_nlcc(i)*upf_pot%r(:)**2)
            END DO
            errcc = sum((cgauss(:) - upf_pot%rho_nlcc(:))**2*upf_pot%r(:)**2*upf_pot%rab(:))
            errcc = sqrt(errcc/real(n, kind=dp))
            DEALLOCATE (cgauss)
         ELSE
            cpabort("sgp_pot%has_nlcc requires is_upf")
         END IF
         error(3) = errcc
      END IF
      !
      IF (is_upf) THEN
         DEALLOCATE (cutpotu)
         DEALLOCATE (cutpots)
      ELSE
         DEALLOCATE (cutpots)
      END IF
      !
      CALL release_opmat(score)
      CALL release_opmat(shnl)
      CALL release_opmat(core)
      CALL release_opmat(hnl)
 
      CALL release_atom_basis(basis)
 
   END SUBROUTINE sgp_construction
 
! **************************************************************************************************
!> \brief ...
!> \param atom_info ...
!> \param input_section ...
!> \param iw ...
! **************************************************************************************************
   SUBROUTINE atom_sgp_construction(atom_info, input_section, iw)
 
      TYPE(atom_p_type), DIMENSION(:, :), POINTER        :: atom_info
      TYPE(section_vals_type), POINTER                   :: input_section
      INTEGER, INTENT(IN)                                :: iw
 
      INTEGER                                            :: i, n, ppot_type
      LOGICAL                                            :: do_transform, explicit, is_ecp, is_upf
      REAL(kind=dp)                                      :: errcc, rcut
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: cgauss, cutpots, cutpotu
      TYPE(atom_ecppot_type), POINTER                    :: ecp_pot
      TYPE(atom_sgp_potential_type)                      :: sgp_pot
      TYPE(atom_type), POINTER                           :: atom_ref
      TYPE(atom_upfpot_type), POINTER                    :: upf_pot
      TYPE(opmat_type), POINTER                          :: core, hnl, score, shnl
 
      CALL section_vals_get(input_section, explicit=explicit)
      IF (.NOT. explicit) RETURN
 
      IF (iw > 0) WRITE (iw, '(/," ",79("*"),/,T24,A,/," ",79("*"))') "SEPARABLE GAUSSIAN PSEUDOPOTENTIAL"
 
      atom_ref => atom_info(1, 1)%atom
 
      ppot_type = atom_ref%potential%ppot_type
      SELECT CASE (ppot_type)
      CASE (gth_pseudo)
         IF (iw > 0) WRITE (iw, '(" GTH Pseudopotential is already in SGP form. ")')
         do_transform = .false.
      CASE (ecp_pseudo)
         do_transform = .true.
         is_ecp = .true.
         is_upf = .false.
         ecp_pot => atom_ref%potential%ecp_pot
      CASE (upf_pseudo)
         do_transform = .true.
         is_ecp = .false.
         is_upf = .true.
         upf_pot => atom_ref%potential%upf_pot
      CASE (no_pseudo)
         IF (iw > 0) WRITE (iw, '(" No Pseudopotential available for transformation. ")')
         do_transform = .false.
      CASE DEFAULT
         cpabort("Unknown ppot_type")
      END SELECT
 
      ! generate the transformed potentials
      IF (do_transform) THEN
         IF (is_ecp) THEN
            CALL ecp_sgp_constr(ecp_pot, sgp_pot, atom_ref%basis)
         ELSE IF (is_upf) THEN
            CALL upf_sgp_constr(upf_pot, sgp_pot, atom_ref%basis)
         ELSE
            cpabort("Either ecp_pseudo or upf_pseudo is needed for atom_sgp_construction")
         END IF
      END IF
 
      ! Check the result
      IF (do_transform) THEN
         NULLIFY (core, hnl)
         CALL create_opmat(core, atom_ref%basis%nbas)
         CALL create_opmat(hnl, atom_ref%basis%nbas, 5)
         NULLIFY (score, shnl)
         CALL create_opmat(score, atom_ref%basis%nbas)
         CALL create_opmat(shnl, atom_ref%basis%nbas, 5)
         !
         ! upf has often very small grids, use a smooth cutoff function
         IF (is_upf) THEN
            n = SIZE(upf_pot%r)
            ALLOCATE (cutpotu(n))
            rcut = maxval(upf_pot%r)
            CALL cutpot(cutpotu, upf_pot%r, rcut, 2.5_dp)
            n = atom_ref%basis%grid%nr
            ALLOCATE (cutpots(n))
            CALL cutpot(cutpots, atom_ref%basis%grid%rad, rcut, 2.5_dp)
         ELSE
            n = atom_ref%basis%grid%nr
            ALLOCATE (cutpots(n))
            cutpots = 1.0_dp
         END IF
         !
         IF (is_ecp) THEN
            CALL ecpints(hnl%op, atom_ref%basis, ecp_pot)
         ELSE IF (is_upf) THEN
            CALL upfints(core%op, hnl%op, atom_ref%basis, upf_pot, cutpotu, sgp_pot%ac_local)
         ELSE
            cpabort("Either ecp_pseudo or upf_pseudo is needed for atom_sgp_construction")
         END IF
         !
         CALL sgpints(score%op, shnl%op, atom_ref%basis, sgp_pot, cutpots)
         !
         IF (sgp_pot%has_local) THEN
            n = min(3, ubound(core%op, 3))
            errcc = maxval(abs(core%op(:, :, 0:n) - score%op(:, :, 0:n)))
            IF (iw > 0) THEN
               WRITE (iw, '(" Local part of pseudopotential")')
               WRITE (iw, '(" Number of basis functions ",T77,i4)') sgp_pot%n_local
               WRITE (iw, '(" Max. abs. error of matrix elements ",T65,f16.8)') errcc
            END IF
         END IF
         IF (sgp_pot%has_nonlocal) THEN
            IF (iw > 0) THEN
               WRITE (iw, '(" Nonlocal part of pseudopotential")')
               WRITE (iw, '(" Max. l-quantum number",T77,i4)') sgp_pot%lmax
               WRITE (iw, '(" Number of projector basis functions ",T77,i4)') sgp_pot%n_nonlocal
               DO i = 0, sgp_pot%lmax
                  errcc = maxval(abs(hnl%op(:, :, i) - shnl%op(:, :, i)))
                  WRITE (iw, '(" Max. abs. error of matrix elements: l=",i2,T69,f12.8)') i, errcc
               END DO
            END IF
         END IF
         IF (sgp_pot%has_nlcc) THEN
            IF (is_upf) THEN
               n = SIZE(upf_pot%r)
               ALLOCATE (cgauss(n))
               cgauss = 0.0_dp
               DO i = 1, sgp_pot%n_nlcc
                  cgauss(:) = cgauss(:) + sgp_pot%c_nlcc(i)*exp(-sgp_pot%a_nlcc(i)*upf_pot%r(:)**2)
               END DO
               errcc = sum((cgauss(:) - upf_pot%rho_nlcc(:))**2*upf_pot%r(:)**2*upf_pot%rab(:))
               errcc = sqrt(errcc/real(n, kind=dp))
               DEALLOCATE (cgauss)
            ELSE
               cpabort("sgp_pot%has_nlcc requires is_upf")
            END IF
            IF (iw > 0) THEN
               WRITE (iw, '(" Non-linear core correction: core density")')
               WRITE (iw, '(" Number of basis functions ",T77,i4)') sgp_pot%n_nlcc
               WRITE (iw, '(" RMS error of core density ",T69,f12.8)') errcc
            END IF
         END IF
         !
         IF (is_upf) THEN
            DEALLOCATE (cutpotu)
            DEALLOCATE (cutpots)
         ELSE
            DEALLOCATE (cutpots)
         END IF
         !
         CALL release_opmat(score)
         CALL release_opmat(shnl)
         CALL release_opmat(core)
         CALL release_opmat(hnl)
      END IF
 
      CALL atom_sgp_release(sgp_pot)
 
      IF (iw > 0) WRITE (iw, '(" ",79("*"))')
 
   END SUBROUTINE atom_sgp_construction
 
! **************************************************************************************************
!> \brief ...
!> \param ecp_pot ...
!> \param sgp_pot ...
!> \param basis ...
! **************************************************************************************************
   SUBROUTINE ecp_sgp_constr(ecp_pot, sgp_pot, basis)
 
      TYPE(atom_ecppot_type)                             :: ecp_pot
      TYPE(atom_sgp_potential_type)                      :: sgp_pot
      TYPE(atom_basis_type)                              :: basis
 
      INTEGER                                            :: i, ia, ir, j, k, l, n, na, nl, nr
      REAL(kind=dp)                                      :: alpha, amx, cmx
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: al, cl, cpot, pgauss
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: cmat, qmat, score, sinv, smat, tmat
      REAL(kind=dp), DIMENSION(:), POINTER               :: rad
 
      sgp_pot%has_nlcc = .false.
      sgp_pot%n_nlcc = 0
      sgp_pot%has_local = .false.
      sgp_pot%n_local = 0
 
      ! transform semilocal potential into a separable local form
      sgp_pot%has_nonlocal = .true.
      !
      ! hardcoded number of projectors (equal for all l values)
      nl = 8
      !
      sgp_pot%n_nonlocal = nl
      sgp_pot%lmax = ecp_pot%lmax
      ALLOCATE (sgp_pot%a_nonlocal(nl))
      ALLOCATE (sgp_pot%h_nonlocal(nl, 0:ecp_pot%lmax))
      ALLOCATE (sgp_pot%c_nonlocal(nl, nl, 0:ecp_pot%lmax))
      !
      ALLOCATE (al(nl), cl(nl))
      ALLOCATE (smat(nl, nl), sinv(nl, nl))
      ALLOCATE (tmat(nl, nl), cmat(nl, nl))
      al = 0.0_dp
      amx = maxval(ecp_pot%bpot)
      cmx = amx/0.15_dp
      cmx = cmx**(1._dp/real(nl - 1, dp))
      cmx = 1._dp/cmx
      DO ir = 1, nl
         al(ir) = amx*cmx**(ir - 1)
      END DO
      !
      sgp_pot%a_nonlocal(1:nl) = al(1:nl)
      !
      nr = basis%grid%nr
      rad => basis%grid%rad
      ALLOCATE (cpot(nr), pgauss(nr))
      DO l = 0, ecp_pot%lmax
         na = basis%nbas(l)
         ALLOCATE (score(na, na), qmat(na, nl))
         cpot = 0._dp
         DO k = 1, ecp_pot%npot(l)
            n = ecp_pot%nrpot(k, l)
            alpha = ecp_pot%bpot(k, l)
            cpot(:) = cpot + ecp_pot%apot(k, l)*rad**(n - 2)*exp(-alpha*rad**2)
         END DO
         DO i = 1, na
            DO j = i, na
               score(i, j) = integrate_grid(cpot, basis%bf(:, i, l), basis%bf(:, j, l), basis%grid)
               score(j, i) = score(i, j)
            END DO
         END DO
         ! overlap basis with projectors
         DO i = 1, nl
            pgauss(:) = exp(-al(i)*rad(:)**2)*rad(:)**l
            DO ia = 1, na
               qmat(ia, i) = integrate_grid(basis%bf(:, ia, l), pgauss(:), basis%grid)
            END DO
         END DO
         qmat = sqrt(fourpi)*qmat
         ! tmat = qmat * score * qmat
         tmat(1:nl, 1:nl) = matmul(transpose(qmat(1:na, 1:nl)), matmul(score(1:na, 1:na), qmat(1:na, 1:nl)))
         smat(1:nl, 1:nl) = matmul(transpose(qmat(1:na, 1:nl)), qmat(1:na, 1:nl))
         CALL get_pseudo_inverse_diag(smat(1:nl, 1:nl), sinv(1:nl, 1:nl), 1.e-10_dp)
         cmat(1:nl, 1:nl) = matmul(sinv(1:nl, 1:nl), matmul(tmat(1:nl, 1:nl), sinv(1:nl, 1:nl)))
         cmat(1:nl, 1:nl) = (cmat(1:nl, 1:nl) + transpose(cmat(1:nl, 1:nl)))*0.5_dp
         !
         ! Residium
         !
         CALL diamat_all(cmat(1:nl, 1:nl), cl(1:nl), .true.)
         !
         sgp_pot%h_nonlocal(1:nl, l) = cl(1:nl)
         sgp_pot%c_nonlocal(1:nl, 1:nl, l) = cmat(1:nl, 1:nl)
         sgp_pot%is_nonlocal(l) = .true.
         !
         DEALLOCATE (score, qmat)
      END DO
      DEALLOCATE (cpot, pgauss)
      DEALLOCATE (al, cl, smat, sinv, tmat, cmat)
 
   END SUBROUTINE ecp_sgp_constr
 
! **************************************************************************************************
!> \brief ...
!> \param upf_pot ...
!> \param sgp_pot ...
!> \param basis ...
! **************************************************************************************************
   SUBROUTINE upf_sgp_constr(upf_pot, sgp_pot, basis)
 
      TYPE(atom_upfpot_type)                             :: upf_pot
      TYPE(atom_sgp_potential_type)                      :: sgp_pot
      TYPE(atom_basis_type)                              :: basis
 
      CHARACTER(len=4)                                   :: ptype
      INTEGER                                            :: ia, ib, ipa, ipb, ir, la, lb, na, nl, &
                                                            np, nr
      LOGICAL                                            :: nl_trans
      REAL(kind=dp)                                      :: cpa, cpb, eee, ei, errloc, rc, x(2), zval
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: al, cl, pgauss, pproa, pprob, vgauss, &
                                                            vloc, ww
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: cmat, qmat, score, sinv, smat, tmat
      TYPE(atom_basis_type)                              :: gbasis
      TYPE(opt_state_type)                               :: ostate
 
      IF (upf_pot%is_ultrasoft .OR. upf_pot%is_paw .OR. upf_pot%is_coulomb) THEN
         sgp_pot%has_nonlocal = .false.
         sgp_pot%n_nonlocal = 0
         sgp_pot%has_local = .false.
         sgp_pot%n_local = 0
         sgp_pot%has_nlcc = .false.
         sgp_pot%n_nlcc = 0
         RETURN
      END IF
 
      ! radial grid
      nr = SIZE(upf_pot%r)
      ! weights for integration
      ALLOCATE (ww(nr))
      ww(:) = upf_pot%r(:)**2*upf_pot%rab(:)
 
      ! start with local potential
      sgp_pot%has_local = .true.
      ! fit local potential to Gaussian form
      ALLOCATE (vloc(nr), vgauss(nr))
      ! smearing of core charge
      zval = upf_pot%zion
      ! Try to find an optimal Gaussian charge distribution
      CALL erffit(sgp_pot%ac_local, upf_pot%vlocal, upf_pot%r, zval)
      sgp_pot%zval = zval
      DO ir = 1, nr
         IF (upf_pot%r(ir) < 1.e-12_dp) THEN
            vgauss(ir) = -sqrt(2.0_dp)*zval/rootpi/sgp_pot%ac_local
         ELSE
            rc = upf_pot%r(ir)/sgp_pot%ac_local/sqrt(2.0_dp)
            vgauss(ir) = -zval/upf_pot%r(ir)*erf(rc)
         END IF
      END DO
      vloc(:) = upf_pot%vlocal(:) - vgauss(:)
      !
      CALL atom_basis_gridrep(basis, gbasis, upf_pot%r, upf_pot%rab)
      !
      nl = 12
      ALLOCATE (al(nl), cl(nl))
      ostate%nf = 0
      ostate%nvar = 2
      x(1) = 1.00_dp !starting point of geometric series
      x(2) = 1.20_dp !factor of series
      ostate%rhoend = 1.e-12_dp
      ostate%rhobeg = 5.e-2_dp
      ostate%maxfun = 1000
      ostate%iprint = 1
      ostate%unit = -1
      ostate%state = 0
      DO
         IF (ostate%state == 2) THEN
            DO ir = 1, nl
               al(ir) = x(1)*x(2)**(ir - 1)
            END DO
            CALL pplocal_error(nl, al, cl, vloc, vgauss, gbasis, upf_pot%r, ww, 1, ostate%f)
         END IF
         IF (ostate%state == -1) EXIT
         CALL powell_optimize(ostate%nvar, x, ostate)
      END DO
      ostate%state = 8
      CALL powell_optimize(ostate%nvar, x, ostate)
      DO ir = 1, nl
         al(ir) = x(1)*x(2)**(ir - 1)
      END DO
      CALL pplocal_error(nl, al, cl, vloc, vgauss, gbasis, upf_pot%r, ww, 1, errloc)
      !
      ALLOCATE (sgp_pot%a_local(nl), sgp_pot%c_local(nl))
      sgp_pot%n_local = nl
      sgp_pot%a_local(1:nl) = al(1:nl)
      sgp_pot%c_local(1:nl) = cl(1:nl)
      DEALLOCATE (vloc, vgauss)
      DEALLOCATE (al, cl)
      CALL release_atom_basis(gbasis)
      !
      ptype = adjustl(trim(upf_pot%pseudo_type))
      IF (ptype(1:2) == "NC" .OR. ptype(1:2) == "US") THEN
         nl_trans = .false.
      ELSE IF (ptype(1:2) == "SL") THEN
         nl_trans = .true.
      ELSE
         cpabort("Pseudopotential type: ["//adjustl(trim(ptype))//"] not known")
      END IF
 
      ! purely local pseudopotentials
      IF (upf_pot%l_max < 0) THEN
         sgp_pot%n_nonlocal = 0
         sgp_pot%lmax = -1
         sgp_pot%has_nonlocal = .false.
      ELSE
         ! Non-local pseudopotential in Gaussian form
         IF (nl_trans) THEN
            sgp_pot%has_nonlocal = .true.
            ! semi local pseudopotential
            ! fit to nonlocal form
            ! get basis representation on UPF grid
            nl = 8
            !
            sgp_pot%n_nonlocal = nl
            sgp_pot%lmax = upf_pot%l_max
            ALLOCATE (sgp_pot%a_nonlocal(nl))
            ALLOCATE (sgp_pot%h_nonlocal(nl, 0:upf_pot%l_max))
            ALLOCATE (sgp_pot%c_nonlocal(nl, nl, 0:upf_pot%l_max))
            !
            ALLOCATE (al(nl), cl(nl))
            ALLOCATE (smat(nl, nl), sinv(nl, nl))
            ALLOCATE (tmat(nl, nl), cmat(nl, nl))
            al = 0.0_dp
            DO ir = 1, nl
               al(ir) = 10.0_dp*0.60_dp**(ir - 1)
            END DO
            !
            sgp_pot%a_nonlocal(1:nl) = al(1:nl)
            !
            CALL atom_basis_gridrep(basis, gbasis, upf_pot%r, upf_pot%rab)
            ALLOCATE (pgauss(nr), vloc(nr))
            DO la = 0, upf_pot%l_max
               IF (la == upf_pot%l_local) cycle
               sgp_pot%is_nonlocal(la) = .true.
               na = gbasis%nbas(la)
               ALLOCATE (score(na, na), qmat(na, nl))
               ! Reference matrix
               vloc(:) = upf_pot%vsemi(:, la + 1) - upf_pot%vlocal(:)
               DO ia = 1, na
                  DO ib = ia, na
                     score(ia, ib) = sum(vloc(:)*gbasis%bf(:, ia, la)*gbasis%bf(:, ib, la)*ww(:))
                     score(ib, ia) = score(ia, ib)
                  END DO
               END DO
               ! overlap basis with projectors
               DO ir = 1, nl
                  pgauss(:) = exp(-al(ir)*upf_pot%r(:)**2)*upf_pot%r(:)**la
                  eee = rootpi/(2._dp**(la + 2)*dfac(2*la + 1))/(2._dp*al(ir))**(la + 1.5_dp)
                  pgauss(:) = pgauss(:)/sqrt(eee)
                  DO ia = 1, na
                     qmat(ia, ir) = sum(gbasis%bf(:, ia, la)*pgauss(:)*ww)
                  END DO
               END DO
               ! tmat = qmat * score * qmat
               tmat(1:nl, 1:nl) = matmul(transpose(qmat(1:na, 1:nl)), matmul(score(1:na, 1:na), qmat(1:na, 1:nl)))
               smat(1:nl, 1:nl) = matmul(transpose(qmat(1:na, 1:nl)), qmat(1:na, 1:nl))
               CALL get_pseudo_inverse_diag(smat(1:nl, 1:nl), sinv(1:nl, 1:nl), 1.e-10_dp)
               cmat(1:nl, 1:nl) = matmul(sinv(1:nl, 1:nl), matmul(tmat(1:nl, 1:nl), sinv(1:nl, 1:nl)))
               cmat(1:nl, 1:nl) = (cmat(1:nl, 1:nl) + transpose(cmat(1:nl, 1:nl)))*0.5_dp
               CALL diamat_all(cmat(1:nl, 1:nl), cl(1:nl), .true.)
               !
               ! get back unnormalized Gaussians
               DO ir = 1, nl
                  ei = rootpi/(2._dp**(la + 2)*dfac(2*la + 1))/(2._dp*al(ir))**(la + 1.5_dp)
                  cmat(ir, 1:nl) = cmat(ir, 1:nl)/sqrt(ei)
               END DO
               sgp_pot%h_nonlocal(1:nl, la) = cl(1:nl)
               sgp_pot%c_nonlocal(1:nl, 1:nl, la) = cmat(1:nl, 1:nl)
               sgp_pot%is_nonlocal(la) = .true.
               DEALLOCATE (score, qmat)
            END DO
            ! SQRT(4PI)
            sgp_pot%c_nonlocal = sgp_pot%c_nonlocal/sqrt(fourpi)
            CALL release_atom_basis(gbasis)
            DEALLOCATE (pgauss, vloc)
            DEALLOCATE (al, cl, smat, sinv, tmat, cmat)
         ELSE
            sgp_pot%has_nonlocal = .true.
            ! non local pseudopotential
            ALLOCATE (pproa(nr), pprob(nr), pgauss(nr))
            np = upf_pot%number_of_proj
            nl = 8
            ALLOCATE (al(nl), cl(nl))
            ALLOCATE (smat(nl, nl), sinv(nl, nl))
            ALLOCATE (tmat(nl, nl), cmat(nl, nl))
            al = 0.0_dp
            cl = 0.0_dp
            DO ir = 1, nl
               al(ir) = 10.0_dp*0.60_dp**(ir - 1)
            END DO
            !
            sgp_pot%lmax = maxval(upf_pot%lbeta(:))
            sgp_pot%n_nonlocal = nl
            ALLOCATE (sgp_pot%a_nonlocal(nl))
            ALLOCATE (sgp_pot%h_nonlocal(nl, 0:sgp_pot%lmax))
            ALLOCATE (sgp_pot%c_nonlocal(nl, nl, 0:sgp_pot%lmax))
            !
            sgp_pot%a_nonlocal(1:nl) = al(1:nl)
            !
            CALL atom_basis_gridrep(basis, gbasis, upf_pot%r, upf_pot%rab)
            DO la = 0, sgp_pot%lmax
               sgp_pot%is_nonlocal(la) = .true.
               na = gbasis%nbas(la)
               ALLOCATE (score(na, na), qmat(na, nl))
               ! Reference matrix
               score = 0.0_dp
               DO ipa = 1, np
                  lb = upf_pot%lbeta(ipa)
                  IF (la /= lb) cycle
                  pproa(:) = upf_pot%beta(:, ipa)
                  DO ipb = 1, np
                     lb = upf_pot%lbeta(ipb)
                     IF (la /= lb) cycle
                     pprob(:) = upf_pot%beta(:, ipb)
                     eee = upf_pot%dion(ipa, ipb)
                     DO ia = 1, na
                        cpa = sum(pproa(:)*gbasis%bf(:, ia, la)*ww(:))
                        DO ib = ia, na
                           cpb = sum(pprob(:)*gbasis%bf(:, ib, la)*ww(:))
                           score(ia, ib) = score(ia, ib) + cpa*eee*cpb
                           score(ib, ia) = score(ia, ib)
                        END DO
                     END DO
                  END DO
               END DO
               ! overlap basis with projectors
               DO ir = 1, nl
                  pgauss(:) = exp(-al(ir)*upf_pot%r(:)**2)*upf_pot%r(:)**la
                  eee = rootpi/(2._dp**(la + 2)*dfac(2*la + 1))/(2._dp*al(ir))**(la + 1.5_dp)
                  pgauss(:) = pgauss(:)/sqrt(eee)
                  DO ia = 1, na
                     qmat(ia, ir) = sum(gbasis%bf(:, ia, la)*pgauss(:)*ww)
                  END DO
               END DO
               ! tmat = qmat * score * qmat
               tmat(1:nl, 1:nl) = matmul(transpose(qmat(1:na, 1:nl)), matmul(score(1:na, 1:na), qmat(1:na, 1:nl)))
               smat(1:nl, 1:nl) = matmul(transpose(qmat(1:na, 1:nl)), qmat(1:na, 1:nl))
               CALL get_pseudo_inverse_diag(smat(1:nl, 1:nl), sinv(1:nl, 1:nl), 1.e-10_dp)
               cmat(1:nl, 1:nl) = matmul(sinv(1:nl, 1:nl), matmul(tmat(1:nl, 1:nl), sinv(1:nl, 1:nl)))
               cmat(1:nl, 1:nl) = (cmat(1:nl, 1:nl) + transpose(cmat(1:nl, 1:nl)))*0.5_dp
               CALL diamat_all(cmat(1:nl, 1:nl), cl(1:nl), .true.)
               !
               ! get back unnormalized Gaussians
               DO ir = 1, nl
                  ei = rootpi/(2._dp**(la + 2)*dfac(2*la + 1))/(2._dp*al(ir))**(la + 1.5_dp)
                  cmat(ir, 1:nl) = cmat(ir, 1:nl)/sqrt(ei)
               END DO
               sgp_pot%h_nonlocal(1:nl, la) = cl(1:nl)
               sgp_pot%c_nonlocal(1:nl, 1:nl, la) = cmat(1:nl, 1:nl)
               sgp_pot%is_nonlocal(la) = .true.
               DEALLOCATE (score, qmat)
            END DO
            ! SQRT(4PI)
            sgp_pot%c_nonlocal = sgp_pot%c_nonlocal/sqrt(fourpi)
            CALL release_atom_basis(gbasis)
            DEALLOCATE (pgauss, pproa, pprob)
            DEALLOCATE (al, cl, smat, sinv, tmat, cmat)
         END IF
      END IF
 
      IF (upf_pot%core_correction) THEN
         sgp_pot%has_nlcc = .true.
      ELSE
         sgp_pot%has_nlcc = .false.
         sgp_pot%n_nlcc = 0
      END IF
 
      ! fit core charge to Gaussian form
      IF (sgp_pot%has_nlcc) THEN
         CALL fit_core_density_gaussians(upf_pot%r, upf_pot%rab, upf_pot%rho_nlcc, al, cl)
         sgp_pot%n_nlcc = SIZE(al)
         ALLOCATE (sgp_pot%a_nlcc, source=al)
         ALLOCATE (sgp_pot%c_nlcc, source=cl)
         DEALLOCATE (al, cl)
      END IF
 
      DEALLOCATE (ww)
 
   END SUBROUTINE upf_sgp_constr
 
! **************************************************************************************************
!> \brief Fit a tabulated core density without forming ill-conditioned normal equations.
!>        Column-normalized, radial least squares use a direct SVD. A broad Gaussian
!>        sequence resolves both the core and its smooth tail. Zero-density samples
!>        beyond the tabulated mesh prevent an unconstrained diffuse continuation.
!> \param r Radial grid in bohr.
!> \param rab Radial mesh Jacobian dr/dj.
!> \param rho True spherical density (without the 4*pi*r**2 prefactor).
!> \param alpha Gaussian exponents in inverse square bohr.
!> \param coeff Coefficients multiplying exp(-alpha*r**2).
! **************************************************************************************************
   SUBROUTINE fit_core_density_gaussians(r, rab, rho, alpha, coeff)
      REAL(dp), INTENT(IN)                               :: r(:), rab(:), rho(:)
      REAL(dp), ALLOCATABLE, INTENT(OUT)                 :: alpha(:), coeff(:)
 
      INTEGER, PARAMETER                                 :: ng = 48, ntail = 200
      REAL(dp), PARAMETER                                :: amax = 300.0_dp, amin = 0.03_dp, &
                                                            rcond = 1.e-12_dp
 
      INTEGER                                            :: i, info, lwork, m, nr, rank
      REAL(dp)                                           :: dr, extent, norms(ng), singular(ng), &
                                                            work_query(1)
      REAL(dp), ALLOCATABLE                              :: matrix(:, :), rhs(:, :), rr(:), sw(:), &
                                                            work(:)
 
      nr = SIZE(r)
      cpassert(nr >= 2 .AND. SIZE(rab) == nr .AND. SIZE(rho) == nr)
      cpassert(all(r >= 0.0_dp) .AND. all(rab > 0.0_dp))
      cpassert(all(r(2:) > r(:nr - 1)))
      m = nr + ntail
      ALLOCATE (alpha(ng), coeff(ng), matrix(m, ng), rhs(m, 1), rr(m), sw(m))
      rr(:nr) = r
      sw(:nr) = r*sqrt(rab)
      ! At 50 bohr even the most diffuse fitting Gaussian is below exp(-75).
      extent = max(r(nr) + 1.0_dp, 50.0_dp)
      dr = (extent - r(nr))/real(ntail, dp)
      DO i = 1, ntail
         rr(nr + i) = r(nr) + real(i, dp)*dr
         sw(nr + i) = rr(nr + i)*sqrt(dr)
      END DO
      sw(m) = sw(m)/sqrt(2.0_dp)
      rhs = 0.0_dp
      rhs(:nr, 1) = rho*sw(:nr)
      DO i = 1, ng
         alpha(i) = amin*exp(log(amax/amin)*real(i - 1, dp)/real(ng - 1, dp))
         matrix(:, i) = exp(-alpha(i)*rr**2)*sw
         norms(i) = sqrt(sum(matrix(:, i)**2))
         cpassert(norms(i) > 0.0_dp)
         matrix(:, i) = matrix(:, i)/norms(i)
      END DO
      CALL dgelss(m, ng, 1, matrix, m, rhs, m, singular, rcond, rank, work_query, -1, info)
      cpassert(info == 0)
      lwork = max(1, ceiling(work_query(1)))
      ALLOCATE (work(lwork))
      CALL dgelss(m, ng, 1, matrix, m, rhs, m, singular, rcond, rank, work, lwork, info)
      IF (info /= 0) cpabort("SVD failed while fitting the UPF core density")
      coeff(:) = rhs(:ng, 1)/norms
   END SUBROUTINE fit_core_density_gaussians
 
! **************************************************************************************************
!> \brief ...
!> \param sgp_pot ...
! **************************************************************************************************
   SUBROUTINE atom_sgp_release(sgp_pot)
 
      TYPE(atom_sgp_potential_type)                      :: sgp_pot
 
      IF (ASSOCIATED(sgp_pot%a_nonlocal)) DEALLOCATE (sgp_pot%a_nonlocal)
      IF (ASSOCIATED(sgp_pot%h_nonlocal)) DEALLOCATE (sgp_pot%h_nonlocal)
      IF (ASSOCIATED(sgp_pot%c_nonlocal)) DEALLOCATE (sgp_pot%c_nonlocal)
 
      IF (ASSOCIATED(sgp_pot%a_local)) DEALLOCATE (sgp_pot%a_local)
      IF (ASSOCIATED(sgp_pot%c_local)) DEALLOCATE (sgp_pot%c_local)
 
      IF (ASSOCIATED(sgp_pot%a_nlcc)) DEALLOCATE (sgp_pot%a_nlcc)
      IF (ASSOCIATED(sgp_pot%c_nlcc)) DEALLOCATE (sgp_pot%c_nlcc)
 
      sgp_pot = atom_sgp_potential_type()
 
   END SUBROUTINE atom_sgp_release
 
! **************************************************************************************************
!> \brief ...
!> \param core ...
!> \param hnl ...
!> \param basis ...
!> \param upf_pot ...
!> \param cutpotu ...
!> \param ac_local ...
! **************************************************************************************************
   SUBROUTINE upfints(core, hnl, basis, upf_pot, cutpotu, ac_local)
      REAL(kind=dp), DIMENSION(:, :, 0:)                 :: core, hnl
      TYPE(atom_basis_type), INTENT(INOUT)               :: basis
      TYPE(atom_upfpot_type)                             :: upf_pot
      REAL(kind=dp), DIMENSION(:), INTENT(IN)            :: cutpotu
      REAL(kind=dp), INTENT(IN)                          :: ac_local
 
      CHARACTER(len=4)                                   :: ptype
      INTEGER                                            :: i, j, k1, k2, la, lb, m, n
      REAL(kind=dp)                                      :: rc, zval
      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, upf_pot%r, upf_pot%rab)
 
      ! local pseudopotential
      core = 0._dp
      n = SIZE(upf_pot%r)
      ALLOCATE (spot(n))
      spot(:) = upf_pot%vlocal(:)
      zval = upf_pot%zion
      DO i = 1, n
         IF (upf_pot%r(i) < 1.e-12_dp) THEN
            spot(i) = spot(i) + sqrt2*zval/rootpi/ac_local
         ELSE
            rc = upf_pot%r(i)/ac_local/sqrt2
            spot(i) = spot(i) + zval/upf_pot%r(i)*erf(rc)
         END IF
      END DO
      spot(:) = spot(:)*cutpotu(:)
 
      CALL numpot_matrix(core, spot, gbasis, 0)
      DEALLOCATE (spot)
 
      hnl = 0._dp
      ptype = adjustl(trim(upf_pot%pseudo_type))
      IF (ptype(1:2) == "NC" .OR. ptype(1:2) == "US") THEN
         ! non local pseudopotential
         n = maxval(gbasis%nbas(:))
         m = upf_pot%number_of_proj
         ALLOCATE (spmat(n, m))
         spmat = 0.0_dp
         DO i = 1, m
            la = upf_pot%lbeta(i)
            DO j = 1, gbasis%nbas(la)
               spmat(j, i) = integrate_grid(upf_pot%beta(:, i), gbasis%bf(:, j, la), gbasis%grid)
            END DO
         END DO
         DO i = 1, m
            la = upf_pot%lbeta(i)
            DO j = 1, m
               lb = upf_pot%lbeta(j)
               IF (la == lb) THEN
                  DO k1 = 1, gbasis%nbas(la)
                     DO k2 = 1, gbasis%nbas(la)
                        hnl(k1, k2, la) = hnl(k1, k2, la) + spmat(k1, i)*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, upf_pot%l_max
            IF (la == upf_pot%l_local) cycle
            m = SIZE(upf_pot%vsemi(:, la + 1))
            ALLOCATE (spot(m))
            spot(:) = upf_pot%vsemi(:, la + 1) - upf_pot%vlocal(:)
            spot(:) = spot(:)*cutpotu(:)
            n = basis%nbas(la)
            DO i = 1, n
               DO j = i, n
                  hnl(i, j, la) = hnl(i, j, la) + &
                                  integrate_grid(spot(:), &
                                                 gbasis%bf(:, i, la), gbasis%bf(:, j, la), gbasis%grid)
                  hnl(j, i, la) = hnl(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 upfints
 
! **************************************************************************************************
!> \brief ...
!> \param hnl ...
!> \param basis ...
!> \param ecp_pot ...
! **************************************************************************************************
   SUBROUTINE ecpints(hnl, basis, ecp_pot)
      REAL(kind=dp), DIMENSION(:, :, 0:)                 :: hnl
      TYPE(atom_basis_type), INTENT(INOUT)               :: basis
      TYPE(atom_ecppot_type)                             :: ecp_pot
 
      INTEGER                                            :: i, j, k, l, m, n
      REAL(kind=dp)                                      :: alpha
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: cpot
      REAL(kind=dp), DIMENSION(:), POINTER               :: rad
 
      rad => basis%grid%rad
      m = basis%grid%nr
      ALLOCATE (cpot(1:m))
 
      ! non local pseudopotential
      hnl = 0.0_dp
      DO l = 0, ecp_pot%lmax
         cpot = 0._dp
         DO k = 1, ecp_pot%npot(l)
            n = ecp_pot%nrpot(k, l)
            alpha = ecp_pot%bpot(k, l)
            cpot(:) = cpot(:) + 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)
               hnl(i, j, l) = integrate_grid(cpot, basis%bf(:, i, l), basis%bf(:, j, l), basis%grid)
               hnl(j, i, l) = hnl(i, j, l)
            END DO
         END DO
      END DO
      DEALLOCATE (cpot)
 
   END SUBROUTINE ecpints
 
! **************************************************************************************************
!> \brief ...
!> \param core ...
!> \param hnl ...
!> \param basis ...
!> \param sgp_pot ...
!> \param cutpots ...
! **************************************************************************************************
   SUBROUTINE sgpints(core, hnl, basis, sgp_pot, cutpots)
      REAL(kind=dp), DIMENSION(:, :, 0:)                 :: core, hnl
      TYPE(atom_basis_type), INTENT(INOUT)               :: basis
      TYPE(atom_sgp_potential_type)                      :: sgp_pot
      REAL(kind=dp), DIMENSION(:), INTENT(IN)            :: cutpots
 
      INTEGER                                            :: i, ia, j, l, m, n, na
      REAL(kind=dp)                                      :: a, 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
      ALLOCATE (cpot(m))
      IF (sgp_pot%has_local) THEN
         zval = sgp_pot%zval
         core = 0._dp
         cpot = 0.0_dp
         DO i = 1, sgp_pot%n_local
            cpot(:) = cpot(:) + sgp_pot%c_local(i)*exp(-sgp_pot%a_local(i)*rad(:)**2)
         END DO
         cpot(:) = cpot(:)*cutpots(:)
         CALL numpot_matrix(core, cpot, basis, 0)
      END IF
      DEALLOCATE (cpot)
 
      ! nonlocal pseudopotential
      IF (sgp_pot%has_nonlocal) THEN
         hnl = 0.0_dp
         ALLOCATE (pgauss(1:m))
         n = sgp_pot%n_nonlocal
         !
         DO l = 0, sgp_pot%lmax
            cpassert(l <= ubound(basis%nbas, 1))
            IF (.NOT. sgp_pot%is_nonlocal(l)) cycle
            ! overlap (a|p)
            na = basis%nbas(l)
            ALLOCATE (qmat(na, n))
            DO i = 1, n
               a = sgp_pot%a_nonlocal(i)
               pgauss(:) = cutpots(:)*exp(-a*rad(:)**2)*rad(:)**l
               DO ia = 1, na
                  qmat(ia, i) = integrate_grid(basis%bf(:, ia, l), pgauss(:), basis%grid)
               END DO
            END DO
            qmat(1:na, 1:n) = sqrt(fourpi)*matmul(qmat(1:na, 1:n), sgp_pot%c_nonlocal(1:n, 1:n, l))
            DO i = 1, na
               DO j = i, na
                  DO ia = 1, n
                     hnl(i, j, l) = hnl(i, j, l) + qmat(i, ia)*qmat(j, ia)*sgp_pot%h_nonlocal(ia, l)
                  END DO
                  hnl(j, i, l) = hnl(i, j, l)
               END DO
            END DO
            DEALLOCATE (qmat)
         END DO
         DEALLOCATE (pgauss)
      END IF
 
   END SUBROUTINE sgpints
 
! **************************************************************************************************
!> \brief ...
!> \param ac ...
!> \param vlocal ...
!> \param r ...
!> \param z ...
! **************************************************************************************************
   SUBROUTINE erffit(ac, vlocal, r, z)
      REAL(kind=dp), INTENT(OUT)                         :: ac
      REAL(kind=dp), DIMENSION(:), INTENT(IN)            :: vlocal, r
      REAL(kind=dp), INTENT(IN)                          :: z
 
      REAL(kind=dp), PARAMETER                           :: rcut = 1.4_dp
 
      INTEGER                                            :: i, j, m, m1
      REAL(kind=dp)                                      :: a1, a2, an, e2, en, rc
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: epot, rval, vpot
 
      m = SIZE(r)
      ALLOCATE (epot(m), vpot(m), rval(m))
      cpassert(SIZE(vlocal) == m)
      IF (r(1) > r(m)) THEN
         DO i = 1, m
            vpot(m - i + 1) = vlocal(i)
            rval(m - i + 1) = r(i)
         END DO
      ELSE
         vpot(1:m) = vlocal(1:m)
         rval(1:m) = r(1:m)
      END IF
      m1 = 1
      DO i = 1, m
         IF (rval(i) > rcut) THEN
            m1 = i
            EXIT
         END IF
      END DO
 
      a1 = 0.2_dp
      a2 = 0.2_dp
      e2 = 1.e20_dp
      epot = 0.0_dp
      DO i = 0, 20
         an = a1 + i*0.025_dp
         rc = 1._dp/(an*sqrt(2.0_dp))
         DO j = m1, m
            epot(j) = vpot(j) + z/rval(j)*erf(rval(j)*rc)
         END DO
         en = sum(abs(epot(m1:m)*rval(m1:m)**2))
         IF (en < e2) THEN
            e2 = en
            a2 = an
         END IF
      END DO
      ac = a2
 
      DEALLOCATE (epot, vpot, rval)
 
   END SUBROUTINE erffit
 
! **************************************************************************************************
!> \brief ...
!> \param nl ...
!> \param al ...
!> \param cl ...
!> \param vloc ...
!> \param vgauss ...
!> \param gbasis ...
!> \param rad ...
!> \param ww ...
!> \param method ...
!> \param errloc ...
! **************************************************************************************************
   SUBROUTINE pplocal_error(nl, al, cl, vloc, vgauss, gbasis, rad, ww, method, errloc)
      INTEGER                                            :: nl
      REAL(kind=dp), DIMENSION(:)                        :: al, cl, vloc, vgauss
      TYPE(atom_basis_type)                              :: gbasis
      REAL(kind=dp), DIMENSION(:)                        :: rad, ww
      INTEGER, INTENT(IN)                                :: method
      REAL(kind=dp)                                      :: errloc
 
      INTEGER                                            :: ia, ib, ir, ix, la, na
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: tv
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: rmat, sinv, smat
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: gmat
 
      cl = 0.0_dp
      IF (method == 1) THEN
         ALLOCATE (tv(nl), smat(nl, nl), sinv(nl, nl))
         DO ir = 1, nl
            vgauss(:) = exp(-al(ir)*rad(:)**2)
            DO ix = 1, nl
               smat(ir, ix) = sum(vgauss(:)*exp(-al(ix)*rad(:)**2)*ww(:))
            END DO
            tv(ir) = sum(vloc(:)*vgauss(:)*ww(:))
         END DO
         CALL get_pseudo_inverse_diag(smat(1:nl, 1:nl), sinv(1:nl, 1:nl), 1.e-12_dp)
         cl(1:nl) = matmul(sinv(1:nl, 1:nl), tv(1:nl))
      ELSE
         !
         ALLOCATE (tv(nl), smat(nl, nl), sinv(nl, nl))
         !
         smat = 0.0_dp
         tv = 0.0_dp
         DO la = 0, min(ubound(gbasis%nbas, 1), 3)
            na = gbasis%nbas(la)
            ALLOCATE (rmat(na, na), gmat(na, na, nl))
            rmat = 0.0_dp
            gmat = 0.0_dp
            DO ia = 1, na
               DO ib = ia, na
                  rmat(ia, ib) = sum(gbasis%bf(:, ia, la)*gbasis%bf(:, ib, la)*vloc(:)*ww(:))
                  rmat(ib, ia) = rmat(ia, ib)
               END DO
            END DO
            DO ir = 1, nl
               vgauss(:) = exp(-al(ir)*rad(:)**2)
               DO ia = 1, na
                  DO ib = ia, na
                     gmat(ia, ib, ir) = sum(gbasis%bf(:, ia, la)*gbasis%bf(:, ib, la)*vgauss(:)*ww(:))
                     gmat(ib, ia, ir) = gmat(ia, ib, ir)
                  END DO
               END DO
            END DO
            DO ir = 1, nl
               tv(ir) = tv(ir) + accurate_dot_product(rmat, gmat(:, :, ir))
               DO ix = ir, nl
                  smat(ir, ix) = smat(ir, ix) + accurate_dot_product(gmat(:, :, ix), gmat(:, :, ir))
                  smat(ix, ir) = smat(ir, ix)
               END DO
            END DO
            DEALLOCATE (rmat, gmat)
         END DO
         sinv = 0.0_dp
         CALL get_pseudo_inverse_diag(smat(1:nl, 1:nl), sinv(1:nl, 1:nl), 1.e-12_dp)
         cl(1:nl) = matmul(sinv(1:nl, 1:nl), tv(1:nl))
      END IF
      !
      vgauss = 0.0_dp
      DO ir = 1, nl
         vgauss(:) = vgauss(:) + cl(ir)*exp(-al(ir)*rad(:)**2)
      END DO
      errloc = sum((vgauss - vloc)**2*ww)
      !
      DEALLOCATE (tv, smat, sinv)
      !
   END SUBROUTINE pplocal_error
 
! **************************************************************************************************
!> \brief ...
!> \param pot ...
!> \param r ...
!> \param rcut ...
!> \param rsmooth ...
! **************************************************************************************************
   SUBROUTINE cutpot(pot, r, rcut, rsmooth)
      REAL(kind=dp), DIMENSION(:), INTENT(INOUT)         :: pot
      REAL(kind=dp), DIMENSION(:), INTENT(IN)            :: r
      REAL(kind=dp), INTENT(IN)                          :: rcut, rsmooth
 
      INTEGER                                            :: i, n
      REAL(kind=dp)                                      :: rab, rx, x
 
      n = SIZE(pot)
      cpassert(n <= SIZE(r))
 
      pot(:) = 1.0_dp
      DO i = 1, n
         rab = r(i)
         IF (rab > rcut) THEN
            pot(i) = 0.0_dp
         ELSE IF (rab > rcut - rsmooth) THEN
            rx = rab - (rcut - rsmooth)
            x = rx/rsmooth
            pot(i) = -6._dp*x**5 + 15._dp*x**4 - 10._dp*x**3 + 1._dp
         END IF
      END DO
 
   END SUBROUTINE cutpot
 
END MODULE atom_sgp