atom_upf_unittest.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                                                      !
!--------------------------------------------------------------------------------------------------!
 
PROGRAM atom_upf_unittest
   USE atom_basis,                      ONLY: atom_basis_opt
   USE atom_kind_orbitals,              ONLY: sgp_potential_conversion
   USE atom_operators,                  ONLY: atom_ppint_release,&
                                              atom_ppint_setup
   USE atom_sgp,                        ONLY: atom_sgp_potential_type,&
                                              atom_sgp_release,&
                                              fit_core_density_gaussians,&
                                              sgp_construction
   USE atom_types,                      ONLY: &
        cgto_basis, atom_basis_gridrep, atom_basis_type, atom_integrals, atom_potential_type, &
        init_atom_basis, init_atom_basis_default_pp, init_atom_potential, release_atom_basis, &
        release_atom_potential, sgp_pseudo, upf_pseudo
   USE atom_upf,                        ONLY: atom_read_upf,&
                                              atom_release_upf,&
                                              atom_upfpot_type
   USE atom_utils,                      ONLY: atom_core_density,&
                                              atom_density
   USE cp_files,                        ONLY: close_file,&
                                              open_file
   USE cp_log_handling,                 ONLY: cp_get_default_logger,&
                                              cp_logger_get_default_io_unit,&
                                              cp_logger_type
   USE cp_output_handling,              ONLY: silent_print_level
   USE external_potential_types,        ONLY: allocate_potential,&
                                              deallocate_potential,&
                                              init_potential,&
                                              sgp_potential_type
   USE ieee_arithmetic,                 ONLY: ieee_is_finite
   USE input_constants,                 ONLY: contracted_gto,&
                                              do_gapw_gct,&
                                              do_gapw_log,&
                                              do_rhf_atom
   USE input_cp2k_atom,                 ONLY: create_atom_section
   USE input_section_types,             ONLY: section_release,&
                                              section_type,&
                                              section_vals_create,&
                                              section_vals_get_subs_vals,&
                                              section_vals_release,&
                                              section_vals_type,&
                                              section_vals_val_set
   USE kinds,                           ONLY: default_path_length,&
                                              dp
   USE libcp2k,                         ONLY: cp2k_finalize,&
                                              cp2k_init
   USE mathconstants,                   ONLY: fourpi,&
                                              pi
   USE orbital_pointers,                ONLY: co,&
                                              init_orbital_pointers,&
                                              nco,&
                                              nso
   USE orbital_transformation_matrices, ONLY: init_spherical_harmonics
   USE qs_grid_atom,                    ONLY: allocate_grid_atom,&
                                              create_grid_atom,&
                                              deallocate_grid_atom,&
                                              grid_atom_type
   USE upf_types,                       ONLY: upf_has_complete_partial_waves
 
   IMPLICIT NONE
   LOGICAL :: source_root_argument
 
   CHARACTER(LEN=default_path_length) :: filename
   INTEGER :: i, io_unit
   REAL(kind=dp) :: error(6)
   TYPE(atom_upfpot_type) :: pot
   TYPE(atom_sgp_potential_type) :: sgp
 
   CALL cp2k_init()
   io_unit = cp_logger_get_default_io_unit()
   CALL check_core_derivatives()
   CALL check_projector_addition_theorem()
   CALL check_atomic_potential()
   CALL check_upf_units()
   CALL check_nlcc_fit()
   CALL check_basis_origin()
   CALL check_basis_derivatives()
   CALL check_radial_grid_endpoints()
   CALL check_contracted_basis_reader()
   CALL check_basis_optimization_export()
   CALL check_uspp_import()
   CALL check_uspp_partial_waves()
 
   ! The same temporary is reused for successive KIND sections in Quickstep.
   ALLOCATE (pot%rho_nlcc(3))
   pot%rho_nlcc = 1.0_dp
   pot%core_correction = .true.
   pot%econf = [2, 4, 0, 0]
   pot%pname = "previous kind"
   CALL atom_release_upf(pot)
   IF (ALLOCATED(pot%rho_nlcc)) error stop "UPF release leaked the core density"
   IF (pot%core_correction) error stop "UPF release retained core-correction metadata"
   IF (any(pot%econf /= 0)) error stop "UPF release retained the previous electron configuration"
   IF (pot%pname /= "") error stop "UPF release retained the previous potential name"
   CALL atom_release_upf(pot)
 
   ALLOCATE (sgp%a_nlcc(2), sgp%c_nlcc(2))
   sgp%has_nlcc = .true.
   sgp%n_nlcc = 2
   sgp%is_nonlocal(2) = .true.
   CALL atom_sgp_release(sgp)
   IF (ASSOCIATED(sgp%a_nlcc) .OR. ASSOCIATED(sgp%c_nlcc)) error stop "SGP release leaked NLCC arrays"
   IF (sgp%has_nlcc .OR. sgp%n_nlcc /= 0) error stop "SGP release retained NLCC metadata"
   IF (any(sgp%is_nonlocal)) error stop "SGP release retained nonlocal channels"
   CALL atom_sgp_release(sgp)
 
   ! Optional external UPF files exercise the real reader and Gaussian transformation.
   ! This also permits auditing a complete SSSP library without bundling it in CP2K.
   DO i = 1, command_argument_count()
      CALL get_command_argument(i, filename)
      ! do_regtest.py passes the CP2K source root before any optional UPF audit files.
      IF (i == 1) THEN
         INQUIRE (file=trim(filename)//"/tests/UNIT_TESTS", exist=source_root_argument)
         IF (source_root_argument) cycle
      END IF
      CALL atom_read_upf(pot, trim(filename))
      IF (abs(sum(pot%econf) - pot%zion) > 1.e-10_dp) error stop "UPF valence count mismatch"
      IF (.NOT. all(ieee_is_finite(pot%vlocal))) error stop "Non-finite UPF local potential"
      IF (any(pot%r(2:) <= pot%r(:SIZE(pot%r) - 1))) error stop "UPF mesh is not increasing"
      IF (any(pot%rab <= 0.0_dp)) error stop "Non-positive UPF mesh weights"
      IF (ALLOCATED(pot%dion)) THEN
         IF (.NOT. all(ieee_is_finite(pot%dion))) error stop "Non-finite UPF projector coefficients"
         IF (.NOT. all(ieee_is_finite(pot%beta))) error stop "Non-finite UPF projectors"
         IF (maxval(abs(pot%dion - transpose(pot%dion))) > 1.e-10_dp) error stop "Non-Hermitian UPF D matrix"
      END IF
      CALL sgp_construction(sgp, upf_pot=pot, error=error)
      IF (.NOT. all(ieee_is_finite(error))) error stop "Non-finite UPF transformation errors"
      IF (pot%core_correction) THEN
         IF (.NOT. ASSOCIATED(sgp%a_nlcc)) error stop "Missing transformed UPF core density"
         IF (.NOT. all(ieee_is_finite(sgp%c_nlcc))) error stop "Non-finite transformed UPF core density"
      END IF
      IF (io_unit > 0) WRITE (io_unit, '(A,1X,A,1X,6ES16.7)') "UPF transformation:", trim(filename), error
      CALL atom_sgp_release(sgp)
      CALL atom_release_upf(pot)
   END DO
 
   CALL cp2k_finalize()
 
CONTAINS
 
! **************************************************************************************************
!> \brief ...
! **************************************************************************************************
   SUBROUTINE check_uspp_import()
      CHARACTER(LEN=*), PARAMETER :: fixture = "uspp-augmentation-test.upf"
 
      INTEGER                                            :: unit
      TYPE(atom_upfpot_type)                             :: pot
      TYPE(cp_logger_type), POINTER                      :: logger
 
      logger => cp_get_default_logger()
      IF (logger%para_env%is_source()) THEN
         OPEN (newunit=unit, file=fixture, status='REPLACE', action='WRITE')
         WRITE (unit, '(A)') &
            '<UPF version="2.0.1">', &
            '<PP_HEADER element="Pt" pseudo_type="USPP" is_ultrasoft="true"', &
            ' z_valence="16" mesh_size="3" number_of_proj="2" l_max="1" />', &
            '<PP_MESH mesh="3">', &
            '<PP_R>', &
            '0.5 1.0 1.5', &
            '</PP_R>', &
            '<PP_RAB>', &
            '0.5 0.5 0.5', &
            '</PP_RAB>', &
            '</PP_MESH>', &
            '<PP_LOCAL type="real" size="3">', &
            '6 4 2', &
            '</PP_LOCAL>', &
            '<PP_NONLOCAL>', &
            '<PP_BETA.1 type="real" size="3" index="1" angular_momentum="0">', &
            '1 4 9', &
            '</PP_BETA.1>', &
            '<PP_BETA.2 type="real" size="3" index="2" angular_momentum="1">', &
            '2 -1 3', &
            '</PP_BETA.2>', &
            '<PP_DIJ>', &
            '2 0 0 8', &
            '</PP_DIJ>', &
            '<PP_AUGMENTATION q_with_l="T" nqf="0" nqlc="3">', &
            '<PP_Q>', &
            '2 0 0 3', &
            '</PP_Q>', &
            '<PP_QIJL.2.2.0 first_index="2" second_index="2" composite_index="3" angular_momentum="0">', &
            '-3 4 5', &
            '</PP_QIJL.2.2.0>', &
            '<PP_QIJL.1.1.0 type="real" size="3" first_index="1" second_index="1"', &
            ' composite_index="1" angular_momentum="0" is_null="false">', &
            '1 2 3', &
            '</PP_QIJL.1.1.0>', &
            '<PP_QIJL.1.2.1 first_index="1" second_index="2" composite_index="2" angular_momentum="1">', &
            '4 -5 6', &
            '</PP_QIJL.1.2.1>', &
            '<PP_QIJL.2.2.2 first_index="2" second_index="2" composite_index="3"', &
            ' angular_momentum="2" is_null="T" size="0"/>', &
            '</PP_AUGMENTATION>', &
            '</PP_NONLOCAL>', &
            '</UPF>'
         CLOSE (unit)
      END IF
      CALL logger%para_env%sync()
      CALL atom_read_upf(pot, fixture, read_augmentation=.true.)
      IF (any(pot%econf(0:3) /= [1, 6, 9, 0])) error stop 'Wrong 62-electron core subtraction'
      IF (any(shape(pot%qijl) /= [3, 3, 3]) .OR. lbound(pot%qijl, 3) /= 0) error stop 'Wrong augmentation shape'
      IF (maxval(abs(pot%qij - reshape([2._dp, 0._dp, 0._dp, 3._dp], [2, 2]))) > 1.e-14_dp) THEN
         error stop 'PP_Q normalization changed'
      END IF
      IF (maxval(abs(pot%qijl(:, 1, 0) - [1._dp, 2._dp, 3._dp])) > 1.e-14_dp) error stop 'Wrong s-s augmentation'
      IF (maxval(abs(pot%qijl(:, 2, 1) - [4._dp, -5._dp, 6._dp])) > 1.e-14_dp) error stop 'Wrong s-p augmentation'
      IF (maxval(abs(pot%qijl(:, 3, 0) - [-3._dp, 4._dp, 5._dp])) > 1.e-14_dp) error stop 'Wrong p-p augmentation'
      IF (any(pot%qijl(:, :, 2) /= 0._dp)) error stop 'Null augmentation was not zero initialized'
      IF (maxval(abs(pot%beta(:, 2) - [4._dp, -1._dp, 2._dp])) > 1.e-14_dp) error stop 'USPP beta normalization changed'
      IF (abs(pot%dion(2, 2) - 4._dp) > 1.e-14_dp) error stop 'USPP D is not in Hartree'
      CALL atom_release_upf(pot)
      IF (ALLOCATED(pot%qij) .OR. ALLOCATED(pot%qijl)) error stop 'Augmentation arrays were not released'
      CALL logger%para_env%sync()
      IF (logger%para_env%is_source()) THEN
         OPEN (newunit=unit, file=fixture, status='OLD')
         CLOSE (unit, status='DELETE')
      END IF
 
   END SUBROUTINE check_uspp_import
 
! **************************************************************************************************
!> \brief ...
! **************************************************************************************************
   SUBROUTINE check_uspp_partial_waves()
      CHARACTER(LEN=*), PARAMETER :: fixture = 'uspp-partial-waves-test.upf'
 
      CHARACTER(LEN=5)                                   :: has_wfc
      INTEGER                                            :: unit, variant
      TYPE(atom_upfpot_type)                             :: potential
      TYPE(cp_logger_type), POINTER                      :: logger
 
      logger => cp_get_default_logger()
      DO variant = 1, 2
         has_wfc = 'false'
         IF (variant == 1) has_wfc = 'true'
         IF (logger%para_env%is_source()) THEN
            OPEN (newunit=unit, file=fixture, status='REPLACE', action='WRITE')
            WRITE (unit, '(A)') &
               '<UPF version="2.0.1">', &
               '<PP_HEADER element="H" pseudo_type="US" is_ultrasoft="true" is_paw="false"', &
               ' has_wfc="'//trim(has_wfc)//'" z_valence="1" mesh_size="3"', &
               ' number_of_proj="2" number_of_wfc="1" l_max="1"/>', &
               '<PP_MESH mesh="3">', &
               '<PP_R>', '0.5 1 1.5', '</PP_R>', &
               '<PP_RAB>', '0.5 0.5 0.5', '</PP_RAB>', '</PP_MESH>', &
               '<PP_LOCAL>', '6 4 2', '</PP_LOCAL>', &
               '<PP_NONLOCAL>', &
               '<PP_BETA.1 index="1" angular_momentum="0" size="3">', '1 4 9', '</PP_BETA.1>', &
               '<PP_BETA.2 index="2" angular_momentum="1" size="3">', '2 -1 3', '</PP_BETA.2>', &
               '<PP_DIJ>', '2 0 0 8', '</PP_DIJ>', &
               '<PP_AUGMENTATION q_with_l="T" nqf="0" nqlc="3">', &
               '<PP_Q>', '2 0 0 3', '</PP_Q>', &
               '<PP_QIJL.1.1.0 first_index="1" second_index="1" composite_index="1" angular_momentum="0">', &
               '1 2 3', '</PP_QIJL.1.1.0>', &
               '<PP_QIJL.1.2.1 first_index="1" second_index="2" composite_index="2" angular_momentum="1">', &
               '4 -5 6', '</PP_QIJL.1.2.1>', &
               '<PP_QIJL.2.2.0 first_index="2" second_index="2" composite_index="3" angular_momentum="0">', &
               '-3 4 5', '</PP_QIJL.2.2.0>', &
               '<PP_QIJL.2.2.2 first_index="2" second_index="2" composite_index="3" angular_momentum="2">', &
               '0 0 0', '</PP_QIJL.2.2.2>', '</PP_AUGMENTATION>', '</PP_NONLOCAL>', &
               '<PP_PSWFC>', '<PP_CHI.1 l="0">', '999 999 999', '</PP_CHI.1>', '</PP_PSWFC>'
            IF (variant == 1) THEN
               ! Reordered, deliberately unnormalized waves distinguish every index and side.
               WRITE (unit, '(A)') &
                  '<PP_FULL_WFC number_of_wfc="2">', &
                  '<PP_PSWFC.2 type="real" size="3" l="1">', '2 3 -4', '</PP_PSWFC.2>', &
                  '<PP_AEWFC.1 index="1" l="0">', '-3 5 -7', '</PP_AEWFC.1>', &
                  '<PP_AEWFC.2 index="2" l="1">', '5 -7 9', '</PP_AEWFC.2>', &
                  '<PP_PSWFC.1 l="0">', '1 -2 3', '</PP_PSWFC.1>', '</PP_FULL_WFC>'
            END IF
            WRITE (unit, '(A)') '</UPF>'
            CLOSE (unit)
         END IF
         CALL logger%para_env%sync()
         IF (variant == 1) THEN
            CALL atom_read_upf(potential, fixture, read_header=.true.)
            IF (ALLOCATED(potential%partial_waves)) error stop 'USPP full waves imported without augmentation request'
            CALL atom_release_upf(potential)
         END IF
         CALL atom_read_upf(potential, fixture, read_augmentation=.true.)
         IF (potential%is_paw) error stop 'USPP full waves were classified as PAW'
         IF (variant == 1) THEN
            IF (.NOT. ALLOCATED(potential%partial_waves)) error stop 'Missing optional USPP full waves'
            associate(waves => potential%partial_waves)
               IF (any(shape(waves%ae_wfc) /= [3, 2])) error stop 'USPP full-wave count differs from beta count'
               IF (any(waves%ae_wfc /= reshape([-3._dp, 5._dp, -7._dp, 5._dp, -7._dp, 9._dp], [3, 2]))) THEN
                  error stop 'USPP AE full-wave order or original normalization changed'
               END IF
               IF (any(waves%ps_wfc /= reshape([1._dp, -2._dp, 3._dp, 2._dp, 3._dp, -4._dp], [3, 2]))) THEN
                  error stop 'USPP PS full waves changed or were confused with atomic guesses'
               END IF
               IF (any(waves%angular_momenta /= [0, 1])) error stop 'Wrong USPP full-wave angular momenta'
            END associate
            IF (.NOT. upf_has_complete_partial_waves(potential)) error stop 'Complete original full waves not recognized'
            DEALLOCATE (potential%partial_waves%ps_wfc)
            IF (upf_has_complete_partial_waves(potential)) error stop 'Incomplete full-wave pair accepted'
         ELSE
            IF (upf_has_complete_partial_waves(potential)) error stop 'Atomic guess waves accepted as complete full waves'
            IF (ALLOCATED(potential%partial_waves)) error stop 'USPP without full waves allocated partial-wave data'
         END IF
         CALL atom_release_upf(potential)
         IF (ALLOCATED(potential%partial_waves)) error stop 'USPP full-wave data not released'
         CALL logger%para_env%sync()
      END DO
      IF (logger%para_env%is_source()) THEN
         OPEN (newunit=unit, file=fixture, status='OLD')
         CLOSE (unit, status='DELETE')
      END IF
   END SUBROUTINE check_uspp_partial_waves
 
! **************************************************************************************************
!> \brief Check the repaired radial maps at the size that exposes endpoint cancellation.
! **************************************************************************************************
   SUBROUTINE check_radial_grid_endpoints()
      INTEGER, PARAMETER :: n = 6400, types(2) = [do_gapw_gct, do_gapw_log]
      REAL(dp), PARAMETER :: inner(2) = [4.83458313724687745e-18_dp, 6.97482911687147684e-18_dp], &
         outer(2) = [2.06843066219244770e17_dp, 57.5213142095328827_dp]
 
      INTEGER                                            :: i
      REAL(dp)                                           :: error
      TYPE(grid_atom_type), POINTER                      :: grid
 
      NULLIFY (grid)
      DO i = 1, SIZE(types)
         CALL allocate_grid_atom(grid)
         CALL create_grid_atom(grid, n, 1, 2, 0, types(i))
         IF (.NOT. all(ieee_is_finite(grid%rad)) .OR. .NOT. all(ieee_is_finite(grid%wr))) THEN
            error stop 'Nonfinite radial grid'
         END IF
         IF (any(grid%rad <= 0.0_dp) .OR. any(grid%wr <= 0.0_dp)) error stop 'Nonpositive radial grid'
         IF (any(grid%rad(:n - 1) <= grid%rad(2:))) error stop 'Unordered radial points'
         ! Independent analytic moments include the r**2 radial measure.
         error = abs(sum(grid%wr*exp(-grid%rad)) - 2.0_dp)
         error = max(error, abs(sum(grid%wr*exp(-grid%rad)/grid%rad) - 1.0_dp))
         error = max(error, abs(sum(grid%wr*exp(-grid%rad2)) - sqrt(pi)/4.0_dp))
         IF (error > 1.e-10_dp) error stop 'Radial moment error'
         ! Independent high-precision endpoints reject clipping rounded cosines.
         IF (abs(grid%rad(n)/inner(i) - 1.0_dp) > 1.e-13_dp .OR. &
             abs(grid%rad(1)/outer(i) - 1.0_dp) > 1.e-13_dp) error stop 'Radial grid endpoint error'
         CALL deallocate_grid_atom(grid)
      END DO
   END SUBROUTINE check_radial_grid_endpoints
 
! **************************************************************************************************
!> \brief Check primitive and contracted Gaussian derivatives on UPF meshes starting at zero.
! **************************************************************************************************
   SUBROUTINE check_basis_origin()
      INTEGER                                            :: i, l, n, variant
      REAL(KIND=dp)                                      :: expected(3), r(6001), rab(6001), scale
      TYPE(atom_basis_type)                              :: basis, radial
 
      CALL init_atom_basis_default_pp(basis)
      n = maxval(basis%nbas)
      ! A long supplied mesh must not generate an unused singular auxiliary grid.
      r = [(0.002_dp*real(i - 1, dp), i=1, SIZE(r))]
      rab = 0.002_dp
      CALL atom_basis_gridrep(basis, radial, r, rab)
      IF (any(radial%grid%rad /= r) .OR. any(radial%grid%wr /= rab*r**2)) THEN
         error stop "Supplied UPF mesh changed during basis evaluation"
      END IF
      IF (maxval(abs(radial%bf(:, 1, 0) - exp(-basis%am(1, 0)*r**2))) > 1.e-14_dp) THEN
         error stop "Incorrect basis values on long UPF mesh"
      END IF
      CALL release_atom_basis(radial)
      DO variant = 1, 2
         scale = real(variant, dp)
         IF (variant == 2) THEN
            basis%basis_type = cgto_basis
            ALLOCATE (basis%cm(n, n, 0:ubound(basis%am, 2)))
            basis%cm = 0.0_dp
            DO i = 1, n
               basis%cm(i, i, :) = scale
            END DO
         END IF
         CALL atom_basis_gridrep(basis, radial, [0.0_dp, 1.e-16_dp, 0.1_dp], [0.1_dp, 0.1_dp, 0.1_dp])
         IF (.NOT. all(ieee_is_finite(radial%bf)) .OR. .NOT. all(ieee_is_finite(radial%dbf)) .OR. &
             .NOT. all(ieee_is_finite(radial%ddbf))) error stop "Nonfinite Gaussian derivatives at the UPF origin"
         DO l = 0, ubound(basis%am, 2)
            DO i = 1, basis%nbas(l)
               expected = 0.0_dp
               IF (l == 0) expected = [1.0_dp, 0.0_dp, -2.0_dp*basis%am(i, l)]
               IF (l == 1) expected(2) = 1.0_dp
               IF (l == 2) expected(3) = 2.0_dp
               IF (maxval(abs([radial%bf(1, i, l), radial%dbf(1, i, l), radial%ddbf(1, i, l)] - &
                              scale*expected)) > 1.e-12_dp) error stop "Incorrect Gaussian origin limit"
            END DO
         END DO
         CALL release_atom_basis(radial)
      END DO
      CALL release_atom_basis(basis)
   END SUBROUTINE check_basis_origin
 
! **************************************************************************************************
!> \brief Differentiate radial basis values and densities independently at finite radii.
! **************************************************************************************************
   SUBROUTINE check_basis_derivatives()
      REAL(dp), PARAMETER :: radius(4) = [0.08_dp, 0.3_dp, 0.9_dp, 2.1_dp], step = 2.e-4_dp
 
      INTEGER                                            :: i, k, l, n, variant
      REAL(dp)                                           :: derivative(5), errors(3), first, lap(5), &
                                                            r(5), reference, rho(5), second
      REAL(dp), ALLOCATABLE                              :: pmat(:, :, :)
      TYPE(atom_basis_type)                              :: basis, radial
 
      CALL init_atom_basis_default_pp(basis)
      n = maxval(basis%nbas)
      ! Nonunit exponents and signed mixing exercise primitives and genuine contractions.
      basis%nbas = 2
      basis%nprim = 2
      basis%am(1, :) = 0.37_dp
      basis%am(2, :) = 2.3_dp
      ALLOCATE (pmat(n, n, 0:ubound(basis%am, 2)))
      pmat = 0.0_dp
      pmat(1, 1, :) = 0.8_dp
      pmat(2, 2, :) = 0.4_dp
      pmat(1, 2, :) = -0.17_dp
      pmat(2, 1, :) = -0.17_dp
      errors = 0.0_dp
      DO variant = 1, 2
         IF (variant == 2) THEN
            basis%basis_type = cgto_basis
            ALLOCATE (basis%cm(n, n, 0:ubound(basis%am, 2)))
            basis%cm = 0.0_dp
            basis%cm(1, 1, :) = 0.7_dp
            basis%cm(2, 1, :) = -0.2_dp
            basis%cm(1, 2, :) = 0.1_dp
            basis%cm(2, 2, :) = 1.3_dp
         END IF
         DO k = 1, SIZE(radius)
            r = radius(k) + step*[-2.0_dp, -1.0_dp, 0.0_dp, 1.0_dp, 2.0_dp]
            CALL atom_basis_gridrep(basis, radial, r, spread(step, 1, 5))
            DO l = 0, ubound(basis%am, 2)
               DO i = 1, 2
                  first = dot_product(radial%bf(:, i, l), [1.0_dp, -8.0_dp, 0.0_dp, 8.0_dp, -1.0_dp])/(12*step)
                  second = dot_product(radial%bf(:, i, l), [-1.0_dp, 16.0_dp, -30.0_dp, 16.0_dp, -1.0_dp])/(12*step**2)
                  errors(1) = max(errors(1), abs(first - radial%dbf(3, i, l))/(1 + abs(first)))
                  errors(2) = max(errors(2), abs(second - radial%ddbf(3, i, l))/(1 + abs(second)))
               END DO
            END DO
            CALL atom_density(rho, pmat, radial, ubound(basis%am, 2), typ="RHO")
            CALL atom_density(derivative, pmat, radial, ubound(basis%am, 2), typ="DER")
            CALL atom_density(lap, pmat, radial, ubound(basis%am, 2), typ="LAP", rr=r*r)
            first = dot_product(rho, [1.0_dp, -8.0_dp, 0.0_dp, 8.0_dp, -1.0_dp])/(12*step)
            second = dot_product(rho, [-1.0_dp, 16.0_dp, -30.0_dp, 16.0_dp, -1.0_dp])/(12*step**2)
            errors(1) = max(errors(1), abs(first - derivative(3))/(1 + abs(first)))
            reference = second + 2*first/radius(k)
            errors(3) = max(errors(3), abs(reference - lap(3))/(1 + abs(reference)))
            CALL release_atom_basis(radial)
         END DO
      END DO
      IF (io_unit > 0) WRITE (io_unit, '(A,3ES24.15)') "Atomic radial derivative/Laplacian scaled errors: ", errors
      DEALLOCATE (pmat)
      CALL release_atom_basis(basis)
      IF (maxval(errors) > 2.e-7_dp) error stop "Atomic radial finite-difference derivatives disagree"
   END SUBROUTINE check_basis_derivatives
 
! **************************************************************************************************
!> \brief Check fitted densities and their derivatives against an analytic reference.
!>        Evaluation points differ from the input mesh, and neither input mesh
!>        nor reference exponents coincide with the internal fitting sequence.
! **************************************************************************************************
   SUBROUTINE check_nlcc_fit()
      INTEGER, PARAMETER                                 :: nr = 1601
      REAL(dp), PARAMETER :: amplitudes(3) = [0.02_dp, 1.7_dp, 0.4_dp], &
         exponents(3) = [0.15_dp, 2.3_dp, 35.0_dp]
 
      INTEGER                                            :: i, mesh
      REAL(dp) :: density_error, derivative, derivative_error, dx, r(nr), rab(nr), radius, &
         reference, reference_derivative, rho(nr), value
      REAL(dp), ALLOCATABLE                              :: alpha(:), coeff(:)
 
      DO mesh = 1, 2
         dx = log(20.0_dp/1.e-6_dp)/real(nr - 1, dp)
         DO i = 1, nr
            IF (mesh == 1) THEN
               r(i) = 20.0_dp*real(i - 1, dp)/real(nr - 1, dp)
               rab(i) = 20.0_dp/real(nr - 1, dp)
            ELSE
               r(i) = 1.e-6_dp*exp(dx*real(i - 1, dp))
               rab(i) = dx*r(i)
            END IF
            rho(i) = sum(amplitudes*exp(-exponents*r(i)**2))
         END DO
         CALL fit_core_density_gaussians(r, rab, rho, alpha, coeff)
         density_error = 0.0_dp
         derivative_error = 0.0_dp
         DO i = 0, 1200
            radius = real(i, dp)*0.017_dp
            reference = sum(amplitudes*exp(-exponents*radius**2))
            reference_derivative = -2.0_dp*radius*sum(exponents*amplitudes*exp(-exponents*radius**2))
            value = sum(coeff*exp(-alpha*radius**2))
            derivative = -2.0_dp*radius*sum(alpha*coeff*exp(-alpha*radius**2))
            density_error = max(density_error, abs(value - reference))
            derivative_error = max(derivative_error, abs(derivative - reference_derivative))
         END DO
         IF (io_unit > 0) WRITE (io_unit, '(A,I2,2ES16.7)') &
            "Core Gaussian fit, density and derivative:", mesh, density_error, derivative_error
         IF (density_error > 2.e-7_dp) error stop "Gaussian core-density fit is inaccurate"
         IF (derivative_error > 3.e-6_dp) error stop "Gaussian core-density derivative is inaccurate"
         rho = 0.0_dp
         CALL fit_core_density_gaussians(r, rab, rho, alpha, coeff)
         IF (any(coeff /= 0.0_dp)) error stop "Zero core density generated nonzero Gaussian coefficients"
      END DO
   END SUBROUTINE check_nlcc_fit
 
! **************************************************************************************************
!> \brief Read an independent miniature UPF fixture and check its numerical convention.
!>        The values are deliberately simple and are not a physical pseudopotential.
!>        Checking beta separately matters: an inverse rescaling of D leaves the NC
!>        operator unchanged but would require rescaling USPP augmentation data too.
! **************************************************************************************************
   SUBROUTINE check_upf_units()
      CHARACTER(LEN=*), PARAMETER :: fixture = &
         'upf-units-test-with-a-potential-filename-longer-than-'// &
         'the-default-eighty-character-string-buffer.upf'
 
      INTEGER                                            :: unit
      TYPE(atom_potential_type)                          :: atomic
      TYPE(cp_logger_type), POINTER                      :: logger
      TYPE(section_type), POINTER                        :: specification
      TYPE(section_vals_type), POINTER                   :: input, potential_section
 
      logger => cp_get_default_logger()
      IF (logger%para_env%is_source()) THEN
         CALL open_file(fixture, unit_number=unit, file_status="REPLACE", file_action="WRITE")
         WRITE (unit, '(A)') &
            '<UPF version="2.0.1">', &
            '<PP_HEADER element="He" pseudo_type="NC" z_valence="2"', &
            ' mesh_size="3" number_of_proj="2" l_max="0" />', &
            '<PP_MESH mesh="3">', &
            '<PP_R type="real" size="3">', &
            '0.5 1.0 1.5', &
            '</PP_R>', &
            '<PP_RAB type="real" size="3">', &
            '0.5 0.5 0.5', &
            '</PP_RAB>', &
            '</PP_MESH>', &
            '<PP_LOCAL type="real" size="3">', &
            '6.0 4.0 2.0', &
            '</PP_LOCAL>', &
            '<PP_NONLOCAL>', &
            '<PP_BETA.1 type="real" size="3" index="1" angular_momentum="0">', &
            '1.0 4.0 9.0', &
            '</PP_BETA.1>', &
            '<PP_BETA.2 type="real" size="3" index="2" angular_momentum="0">', &
            '2.0 -1.0 3.0', &
            '</PP_BETA.2>', &
            '<PP_DIJ type="real" size="4">', &
            '2.0 0.4 0.4 8.0', &
            '</PP_DIJ>', &
            '</PP_NONLOCAL>', &
            '</UPF>'
         CALL close_file(unit)
      END IF
      CALL logger%para_env%sync()
      NULLIFY (specification, input)
      CALL create_atom_section(specification)
      CALL section_vals_create(input, specification)
      potential_section => section_vals_get_subs_vals(input, 'POTENTIAL')
      CALL section_vals_val_set(potential_section, 'PSEUDO_TYPE', i_val=upf_pseudo)
      CALL section_vals_val_set(potential_section, 'POTENTIAL_FILE_NAME', c_val=fixture)
      CALL init_atom_potential(atomic, potential_section, 2)
      associate(potential => atomic%upf_pot)
         IF (maxval(abs(potential%vlocal - [3.0_dp, 2.0_dp, 1.0_dp])) > 1.e-14_dp) THEN
            error stop "UPF local potential is not in Hartree"
         END IF
         IF (maxval(abs(potential%beta(:, 1) - [2.0_dp, 4.0_dp, 6.0_dp])) > 1.e-14_dp) THEN
            error stop "UPF changed the first tabulated projector amplitude"
         END IF
         IF (maxval(abs(potential%beta(:, 2) - [4.0_dp, -1.0_dp, 2.0_dp])) > 1.e-14_dp) THEN
            error stop "UPF changed the second tabulated projector amplitude"
         END IF
         IF (maxval(abs(potential%dion - reshape([1.0_dp, 0.2_dp, 0.2_dp, 4.0_dp], [2, 2]))) > 1.e-14_dp) THEN
            error stop "UPF D matrix is not in the tabulated-projector convention"
         END IF
      END associate
      CALL release_atom_potential(atomic)
      CALL section_vals_release(input)
      CALL section_release(specification)
      CALL logger%para_env%sync()
      IF (logger%para_env%is_source()) THEN
         CALL open_file(fixture, unit_number=unit, file_status="OLD")
         CALL close_file(unit, file_status="DELETE")
      END IF
   END SUBROUTINE check_upf_units
 
! **************************************************************************************************
!> \brief Compare the converted atomic SGP operator with analytic Gaussian integrals.
!>        This checks that atomic guesses and the Quickstep Hamiltonian describe
!>        the same radial projectors and erf-regularized local Coulomb potential.
! **************************************************************************************************
   SUBROUTINE check_atomic_potential()
      INTEGER                                            :: l
      REAL(KIND=dp)                                      :: a, b, p, reference
      TYPE(atom_basis_type)                              :: basis
      TYPE(atom_integrals)                               :: integrals
      TYPE(atom_potential_type)                          :: atomic
      TYPE(sgp_potential_type), POINTER                  :: potential
 
      NULLIFY (potential)
      CALL allocate_potential(potential)
      potential%zeff = 1.0_dp
      potential%core_charge_radius = 0.4_dp
      potential%alpha_core_charge = 0.5_dp/potential%core_charge_radius**2
      potential%n_nonlocal = 1
      potential%lmax = 2
      potential%is_nonlocal(0:2) = .true.
      potential%n_local = 1
      ALLOCATE (potential%elec_conf(0:3), potential%a_nonlocal(1), &
                potential%c_nonlocal(1, 1, 0:2), potential%h_nonlocal(1, 0:2), &
                potential%a_local(1), potential%c_local(1))
      potential%elec_conf = [1, 0, 0, 0]
      potential%a_nonlocal = 0.7_dp
      potential%c_nonlocal = 1.0_dp/sqrt(fourpi)
      potential%h_nonlocal = 1.0_dp
      potential%a_local = 1.2_dp
      potential%c_local = 0.3_dp
      atomic%ppot_type = sgp_pseudo
      CALL sgp_potential_conversion(potential, atomic%sgp_pot)
      CALL init_atom_basis_default_pp(basis)
      CALL atom_ppint_setup(integrals, basis, atomic)
      DO l = 0, 2
         a = basis%am(8, l) + potential%a_nonlocal(1)
         reference = (gamma(l + 1.5_dp)/(2.0_dp*a**(l + 1.5_dp)))**2
         IF (abs(integrals%hnl(8, 8, l) - reference) > 1.e-9_dp*abs(reference)) THEN
            error stop "Atomic SGP nonlocal operator differs from analytic Gaussian integrals"
         END IF
      END DO
      p = 2.0_dp*basis%am(8, 0)
      b = 1.0_dp/(sqrt(2.0_dp)*potential%core_charge_radius)
      reference = -potential%zeff*b/(2.0_dp*p*sqrt(p + b*b)) + &
                  potential%c_local(1)*gamma(1.5_dp)/(2.0_dp*(p + potential%a_local(1))**1.5_dp)
      IF (abs(integrals%core(8, 8, 0) - reference) > 1.e-9_dp*abs(reference)) THEN
         error stop "Atomic SGP local operator differs from analytic Gaussian integrals"
      END IF
      CALL atom_ppint_release(integrals)
      CALL release_atom_basis(basis)
      CALL deallocate_potential(potential)
   END SUBROUTINE check_atomic_potential
 
! **************************************************************************************************
!> \brief Verify the physical angular kernel of Gaussian UPF projectors through l=5.
!>        For beta_l(r) = r**l exp(-a*r**2), stripping off the radial exponentials
!>        at unit vectors leaves sum_m Y_lm(u) Y_lm(v) = (2*l+1) P_l(u.v)/(4*pi).
!>        The reference uses the Legendre recurrence, independently of the Cartesian
!>        to spherical transformation used by the production projector integrals.
! **************************************************************************************************
   SUBROUTINE check_projector_addition_theorem()
      INTEGER, PARAMETER                                 :: lmax = 5
 
      INTEGER                                            :: i, j, l, lx, ly, lz, nc, ns
      REAL(KIND=dp)                                      :: cross, diagonal, p(0:lmax), u(3), v(3), x
      REAL(KIND=dp), ALLOCATABLE                         :: pu(:), pv(:)
      TYPE(sgp_potential_type), POINTER                  :: potential
 
      NULLIFY (potential)
      CALL init_orbital_pointers(lmax)
      CALL init_spherical_harmonics(lmax, -1)
      CALL allocate_potential(potential)
      potential%n_nonlocal = 1
      potential%lmax = lmax
      potential%nppnl = sum(nco(0:lmax))
      ALLOCATE (potential%c_nonlocal(1, 1, 0:lmax), potential%h_nonlocal(1, 0:lmax))
      ! atom_sgp stores the radial coefficients divided by sqrt(4*pi).
      potential%c_nonlocal = 1.0_dp/sqrt(fourpi)
      potential%h_nonlocal = 1.0_dp
      CALL init_potential(potential)
      u = [0.6_dp, 0.0_dp, 0.8_dp]
      v = [0.3_dp, -0.4_dp, 0.5_dp]
      v = v/norm2(v)
      x = dot_product(u, v)
      p(0) = 1.0_dp
      p(1) = x
      DO l = 2, lmax
         p(l) = ((2*l - 1)*x*p(l - 1) - (l - 1)*p(l - 2))/l
      END DO
      nc = 0
      ns = 0
      DO l = 0, lmax
         ALLOCATE (pu(nso(l)), pv(nso(l)))
         pu = 0.0_dp
         pv = 0.0_dp
         DO lx = 0, l
            DO ly = 0, l - lx
               lz = l - lx - ly
               i = nc + co(lx, ly, lz)
               DO j = 1, nso(l)
                  pu(j) = pu(j) + potential%cprj_ppnl(i, ns + j)*u(1)**lx*u(2)**ly*u(3)**lz
                  pv(j) = pv(j) + potential%cprj_ppnl(i, ns + j)*v(1)**lx*v(2)**ly*v(3)**lz
               END DO
            END DO
         END DO
         diagonal = real(2*l + 1, dp)/fourpi
         cross = diagonal*p(l)
         IF (abs(dot_product(pu, pu) - diagonal) > 1.e-12_dp .OR. &
             abs(dot_product(pv, pv) - diagonal) > 1.e-12_dp .OR. &
             abs(dot_product(pu, pv) - cross) > 1.e-12_dp) THEN
            error stop "SGP projector kernel violates the spherical-harmonic addition theorem"
         END IF
         nc = nc + nco(l)
         ns = ns + nso(l)
         DEALLOCATE (pu, pv)
      END DO
      CALL deallocate_potential(potential)
   END SUBROUTINE check_projector_addition_theorem
 
! **************************************************************************************************
!> \brief Verify radial derivatives and the spherical Laplacian against finite differences.
! **************************************************************************************************
   SUBROUTINE check_core_derivatives()
      REAL(kind=dp), PARAMETER                           :: h = 1.e-4_dp
 
      REAL(kind=dp)                                      :: deriv(4), lap(4), lap_fd(4), minus(4), &
                                                            plus(4), r(4), rho(4)
      TYPE(atom_potential_type)                          :: potential
 
      potential%ppot_type = sgp_pseudo
      potential%sgp_pot%has_nlcc = .true.
      potential%sgp_pot%n_nlcc = 2
      potential%sgp_pot%a_nlcc(1:2) = [0.7_dp, 2.3_dp]
      potential%sgp_pot%c_nlcc(1:2) = [1.8_dp, -0.2_dp]
      r = [0.0_dp, 0.1_dp, 0.7_dp, 2.0_dp]
      rho = 0.0_dp
      plus = 0.0_dp
      minus = 0.0_dp
      deriv = 0.0_dp
      lap = 0.0_dp
      CALL atom_core_density(rho, potential, "RHO", r)
      CALL atom_core_density(plus, potential, "RHO", r + h)
      CALL atom_core_density(minus, potential, "RHO", r - h)
      CALL atom_core_density(deriv, potential, "DER", r)
      CALL atom_core_density(lap, potential, "LAP", r)
      IF (maxval(abs(deriv - (plus - minus)/(2.0_dp*h))) > 1.e-7_dp) THEN
         error stop "SGP core-density derivative differs from finite differences"
      END IF
      lap_fd = (plus - 2.0_dp*rho + minus)/h**2
      lap_fd(1) = 3.0_dp*lap_fd(1)
      lap_fd(2:) = lap_fd(2:) + (plus(2:) - minus(2:))/(h*r(2:))
      IF (maxval(abs(lap - lap_fd)) > 1.e-6_dp) error stop "SGP core-density Laplacian differs from finite differences"
   END SUBROUTINE check_core_derivatives
 
! **************************************************************************************************
!> \brief Check independent shell, primitive and set counts in file and inline atomic bases.
! **************************************************************************************************
   SUBROUTINE check_contracted_basis_reader()
      INTEGER, PARAMETER :: contractions(3) = [11, 1, 1], lower_l(3) = [0, 1, 2], &
         primitives(3) = [11, 41, 2], sets(3) = [1, 1, 21], upper_l(3) = [3, 6, 3]
      CHARACTER(LEN=1600), ALLOCATABLE                   :: lines(:)
      INTEGER                                            :: b, c, i, j, l, mode, n, nc, np, ns, row, &
                                                            s, t, unit
      REAL(dp)                                           :: alpha, coefficient, expected, &
                                                            normalization
      TYPE(atom_basis_type)                              :: basis
      CHARACTER(LEN=*), PARAMETER :: fixture = "atom-contracted-basis-test"
 
      TYPE(cp_logger_type), POINTER                      :: logger
      TYPE(section_type), POINTER                        :: specification
      TYPE(section_vals_type), POINTER                   :: input, pp_basis
 
      logger => cp_get_default_logger()
      DO b = 1, SIZE(sets)
         np = primitives(b)
         nc = contractions(b)
         ns = sets(b)
         ALLOCATE (lines(1 + ns*(np + 1)))
         WRITE (lines(1), '(I0)') ns
         row = 1
         DO s = 1, ns
            row = row + 1
            WRITE (lines(row), '(*(I0,1X))') 1, lower_l(b), upper_l(b), np, (nc, l=lower_l(b), upper_l(b))
            DO i = 1, np
               row = row + 1
               alpha = real(10*s + i, dp)/20.0_dp
               WRITE (lines(row), '(*(F12.6,1X))') alpha, &
                  ((real(100*s + 10*i + c + 5*l, dp)/1000.0_dp, c=1, nc), l=lower_l(b), upper_l(b))
            END DO
         END DO
         IF (logger%para_env%is_source()) THEN
            OPEN (newunit=unit, file=fixture, status="REPLACE", action="WRITE")
            WRITE (unit, '(A)') 'H READER-TEST', (trim(lines(i)), i=1, SIZE(lines))
            CLOSE (unit)
         END IF
         CALL logger%para_env%sync()
         DO mode = 1, 2
            NULLIFY (specification, input)
            CALL create_atom_section(specification)
            CALL section_vals_create(input, specification)
            pp_basis => section_vals_get_subs_vals(input, 'PP_BASIS')
            CALL section_vals_val_set(pp_basis, 'BASIS_TYPE', i_val=contracted_gto)
            CALL section_vals_val_set(pp_basis, 'GRID_POINTS', i_val=32)
            CALL section_vals_val_set(pp_basis, 'BASIS_SET', c_val='READER-TEST')
            CALL section_vals_val_set(pp_basis, 'BASIS_SET_FILE_NAME', c_val=fixture)
            IF (mode == 2) THEN
               DO i = 1, SIZE(lines)
                  CALL section_vals_val_set(pp_basis, 'BASIS%_DEFAULT_KEYWORD_', i_rep_val=i, c_val=trim(lines(i)))
               END DO
            END IF
            CALL init_atom_basis(basis, pp_basis, 1, 'PP')
            DO l = 0, ubound(basis%nbas, 1)
               n = 0
               IF (l >= lower_l(b) .AND. l <= upper_l(b)) n = ns
               IF (basis%nprim(l) /= n*np .OR. basis%nbas(l) /= n*nc) THEN
                  error stop "Incorrect atomic contracted basis dimensions"
               END IF
               DO s = 1, n
                  DO i = 1, np
                     alpha = real(10*s + i, dp)/20.0_dp
                     j = (s - 1)*np + i
                     IF (abs(basis%am(j, l) - alpha) > 1.e-14_dp) error stop "Incorrect contracted basis exponent"
                     ! Inverse square root of the analytic radial Gaussian norm.
                     normalization = sqrt(2.0_dp*(2.0_dp*alpha)**(l + 1.5_dp)/gamma(l + 1.5_dp))
                     DO t = 1, n
                        DO c = 1, nc
                           coefficient = 0.0_dp
                           IF (s == t) coefficient = real(100*s + 10*i + c + 5*l, dp)/1000.0_dp
                           expected = normalization*coefficient
                           IF (abs(basis%cm(j, (t - 1)*nc + c, l) - expected) > 1.e-12_dp*max(1.0_dp, abs(expected))) THEN
                              error stop "Incorrect contracted basis coefficient or normalization"
                           END IF
                        END DO
                     END DO
                  END DO
               END DO
            END DO
            CALL release_atom_basis(basis)
            CALL section_vals_release(input)
            CALL section_release(specification)
         END DO
         CALL logger%para_env%sync()
         IF (logger%para_env%is_source()) THEN
            OPEN (newunit=unit, file=fixture, status="OLD")
            CLOSE (unit, status="DELETE")
         END IF
         DEALLOCATE (lines)
      END DO
   END SUBROUTINE check_contracted_basis_reader
 
! **************************************************************************************************
!> \brief Verify exported orbitals use the best basis restored after Powell optimization.
! **************************************************************************************************
   SUBROUTINE check_basis_optimization_export()
      CHARACTER(LEN=160)                                 :: line
      INTEGER                                            :: i, j, status, unit
      REAL(dp)                                           :: alpha(2), coefficient(2), norm
      REAL(dp), POINTER                                  :: exponents(:)
      TYPE(cp_logger_type), POINTER                      :: logger
      TYPE(section_type), POINTER                        :: specification
      TYPE(section_vals_type), POINTER                   :: input
 
      NULLIFY (specification, input)
      CALL create_atom_section(specification)
      CALL section_vals_create(input, specification)
      CALL section_vals_val_set(input, 'ELEMENT', c_val='He')
      CALL section_vals_val_set(input, 'ELECTRON_CONFIGURATION', c_val='1s2')
      CALL section_vals_val_set(input, 'MAX_ANGULAR_MOMENTUM', i_val=0)
      CALL section_vals_val_set(input, 'AE_BASIS%NUM_GTO', i_val=2)
      ALLOCATE (exponents(2))
      exponents = [2.0_dp, 0.4_dp]
      CALL section_vals_val_set(input, 'AE_BASIS%S_EXPONENTS', r_vals_ptr=exponents)
      CALL section_vals_val_set(input, 'METHOD%METHOD_TYPE', i_val=do_rhf_atom)
      CALL section_vals_val_set(input, 'OPTIMIZATION%EPS_SCF', r_val=1.e-10_dp)
      CALL section_vals_val_set(input, 'OPTIMIZATION%N_DIIS', i_val=8)
      CALL section_vals_val_set(input, 'OPTIMIZATION%MAX_ITER', i_val=100)
      CALL section_vals_val_set(input, 'POWELL%MAX_FUN', i_val=5)
      CALL section_vals_val_set(input, 'POWELL%STEP_SIZE', r_val=0.15_dp)
      CALL section_vals_val_set(input, 'POWELL%ACCURACY', r_val=1.e-5_dp)
      CALL section_vals_val_set(input, 'PRINT%FIT_BASIS%_SECTION_PARAMETERS_', i_val=silent_print_level)
      CALL atom_basis_opt(input)
      CALL section_vals_release(input)
      CALL section_release(specification)
 
      logger => cp_get_default_logger()
      CALL logger%para_env%sync()
      OPEN (newunit=unit, file='OPT_BASIS', status='OLD', action='READ')
      DO
         READ (unit, '(A)', iostat=status) line
         IF (status /= 0) error stop 'Missing optimized basis orbital coefficients'
         IF (index(line, 'L Quantum Number:') > 0) EXIT
      END DO
      DO i = 1, 2
         READ (unit, *) alpha(i), coefficient(i)
      END DO
      CLOSE (unit)
      IF (alpha(1) <= 2.0_dp) error stop 'Optimization did not restore the earlier best exponent'
      norm = 0.0_dp
      DO i = 1, 2
         DO j = 1, 2
            norm = norm + coefficient(i)*coefficient(j)* &
                   (2.0_dp*sqrt(alpha(i)*alpha(j))/(alpha(i) + alpha(j)))**1.5_dp
         END DO
      END DO
      ! The export rounds its coefficients to eight decimal places.
      IF (.NOT. ieee_is_finite(norm) .OR. abs(norm - 1.0_dp) > 1.e-7_dp) THEN
         error stop 'Exported orbital is inconsistent with the optimized basis'
      END IF
      CALL logger%para_env%sync()
      IF (logger%para_env%is_source()) THEN
         OPEN (newunit=unit, file='OPT_BASIS', status='OLD')
         CLOSE (unit, status='DELETE')
      END IF
   END SUBROUTINE check_basis_optimization_export
 
END PROGRAM atom_upf_unittest