gw_optimize_ri_rs_grid.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief Local-environment optimizer for atom-centred real-space RI grids.
!> \author Jan Wilhelm
! **************************************************************************************************
MODULE gw_optimize_ri_rs_grid
   USE basis_set_types,                 ONLY: gto_basis_set_type
   USE cell_types,                      ONLY: cell_type,&
                                              get_cell,&
                                              pbc,&
                                              use_perd_none
   USE cp_lbfgs,                        ONLY: setulb
   USE gw_integrals,                    ONLY: build_3c_integral_block_ctx,&
                                              gw_3c_ctx_create,&
                                              gw_3c_ctx_release,&
                                              gw_3c_ctx_type,&
                                              gw_3c_ws_create,&
                                              gw_3c_ws_release,&
                                              gw_3c_ws_type
   USE input_constants,                 ONLY: do_potential_truncated
   USE kinds,                           ONLY: dp
   USE libint_2c_3c,                    ONLY: libint_potential_type
   USE message_passing,                 ONLY: mp_para_env_type
   USE particle_types,                  ONLY: particle_type
   USE physcon,                         ONLY: angstrom
   USE post_scf_bandstructure_types,    ONLY: post_scf_bandstructure_type,&
                                              rirs_grid_type
   USE qs_environment_types,            ONLY: qs_environment_type
   USE rirs_grid_utils,                 ONLY: evaluate_ao_basis_on_points
#include "./base/base_uses.f90"
 
   IMPLICIT NONE
   PRIVATE
 
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gw_optimize_ri_rs_grid'
   REAL(KIND=dp), PARAMETER, PRIVATE :: optimizer_accuracy = 1.0e-9_dp
   REAL(KIND=dp), PARAMETER, PRIVATE :: truncated_coulomb_cutoff = 3.0_dp/angstrom
 
   TYPE :: local_cluster_type
      INTEGER, ALLOCATABLE                 :: atoms(:)
      REAL(KIND=dp), ALLOCATABLE           :: three_center(:, :, :)
   END TYPE local_cluster_type
 
   PUBLIC :: optimize_ri_rs_grid
 
CONTAINS
 
! **************************************************************************************************
!> \brief Optimize the selected RI-RS grids stored in a GW band-structure environment.
!> \param qs_env Quickstep environment used by the upstream three-center-integral context.
!> \param bs_env GW environment containing all configuration, basis, geometry, parallel, and grid data.
! **************************************************************************************************
   SUBROUTINE optimize_ri_rs_grid(qs_env, bs_env)
      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
 
      CHARACTER(len=*), PARAMETER :: routinen = 'optimize_ri_rs_grid'
 
      INTEGER                                            :: handle, iatom, ikind, n_variable, &
                                                            periodic(3), unit_nr
      INTEGER, ALLOCATABLE                               :: ao_size(:), grid_offsets(:), ri_size(:)
      LOGICAL                                            :: successful
      REAL(kind=dp)                                      :: f, maximum_absolute_error
      REAL(kind=dp), ALLOCATABLE                         :: g(:), lower(:), upper(:), x(:)
      TYPE(cell_type), POINTER                           :: cell
      TYPE(gw_3c_ctx_type)                               :: integral_context
      TYPE(local_cluster_type), ALLOCATABLE              :: clusters(:)
      TYPE(mp_para_env_type), POINTER                    :: para_env
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
 
      CALL timeset(routinen, handle)
      NULLIFY (cell, para_env, particle_set)
      cell => bs_env%ri_rs%cell
      particle_set => bs_env%ri_rs%particle_set
      para_env => bs_env%para_env
      cpassert(ASSOCIATED(cell))
      cpassert(ASSOCIATED(particle_set))
      cpassert(ASSOCIATED(para_env))
      cpassert(ALLOCATED(bs_env%basis_set_AO))
      cpassert(ALLOCATED(bs_env%basis_set_RI))
      cpassert(ALLOCATED(bs_env%ri_rs%grid_cache))
 
      IF (bs_env%ri_rs%grid_opt%max_iter < 1) THEN
         cpabort("GRID_OPTIMIZATION%MAX_ITER must be positive")
      END IF
      IF (bs_env%ri_rs%grid_opt%cutoff_atomic_cluster <= 0.0_dp) THEN
         cpabort("GRID_OPTIMIZATION%CUTOFF_ATOMIC_CLUSTER must be positive")
      END IF
      IF (bs_env%ri_rs%tikhonov < 0.0_dp) THEN
         cpabort("RI_RS%TIKHONOV must not be negative")
      END IF
 
      CALL get_cell(cell, periodic=periodic)
      IF (any(periodic /= use_perd_none)) THEN
         cpabort("GRID_OPTIMIZATION currently supports nonperiodic local environments only")
      END IF
 
      unit_nr = bs_env%unit_nr
      ALLOCATE (ao_size(SIZE(particle_set)), ri_size(SIZE(particle_set)))
      DO iatom = 1, SIZE(particle_set)
         ikind = particle_set(iatom)%atomic_kind%kind_number
         cpassert(ASSOCIATED(bs_env%basis_set_AO(ikind)%gto_basis_set))
         cpassert(ASSOCIATED(bs_env%basis_set_RI(ikind)%gto_basis_set))
         ao_size(iatom) = bs_env%basis_set_AO(ikind)%gto_basis_set%nsgf
         ri_size(iatom) = bs_env%basis_set_RI(ikind)%gto_basis_set%nsgf
      END DO
      IF (any(ao_size < 1) .OR. any(ri_size < 1)) THEN
         cpabort("Every atom in GRID_OPTIMIZATION needs ORB and RI_AUX functions")
      END IF
 
      cpassert(SIZE(bs_env%ri_rs%grid_cache) == SIZE(particle_set))
      ALLOCATE (grid_offsets(SIZE(bs_env%ri_rs%grid_cache)))
      n_variable = 0
      DO iatom = 1, SIZE(bs_env%ri_rs%grid_cache)
         grid_offsets(iatom) = n_variable
         n_variable = n_variable + 3*SIZE(bs_env%ri_rs%grid_cache(iatom)%raw_points, 2)
      END DO
      block
         TYPE(libint_potential_type) :: potential
         potential%potential_type = do_potential_truncated
         potential%cutoff_radius = truncated_coulomb_cutoff
         potential%omega = 0.0_dp
         potential%filename = "t_c_g.dat"
         CALL gw_3c_ctx_create(integral_context, qs_env, potential, &
                               bs_env%basis_set_AO, bs_env%basis_set_AO, bs_env%basis_set_RI)
      END block
      CALL build_local_clusters(integral_context, particle_set, cell, ao_size, ri_size, &
                                bs_env%ri_rs%grid_opt%cutoff_atomic_cluster, para_env, clusters)
 
      ALLOCATE (x(n_variable), g(n_variable), lower(n_variable), upper(n_variable))
      CALL pack_atom_grids(bs_env%ri_rs%grid_cache, grid_offsets, x)
      lower(:) = x - 100.0_dp/angstrom
      upper(:) = x + 100.0_dp/angstrom
      CALL optimize_grid_coordinates(x, lower, upper, optimizer_accuracy, &
                                     bs_env%ri_rs%grid_opt%max_iter, bs_env, clusters, grid_offsets)
      CALL grid_objective(x, bs_env, clusters, grid_offsets, f, g, &
                          maximum_absolute_error, successful)
      IF (.NOT. successful) cpabort("RI-RS grid optimization produced no regularized fit")
 
      IF (unit_nr > 0) THEN
         WRITE (unit_nr, '(/,T2,A)') 'RI-RS grid optimization terminated'
         WRITE (unit_nr, '(T2,A,T72,ES9.1)') 'Normalized 3C error:', f
         WRITE (unit_nr, '(T2,A,T72,ES9.1)') &
            'Maximum absolute 3C error:', maximum_absolute_error
      END IF
 
      DO iatom = 1, SIZE(bs_env%ri_rs%grid_cache)
         bs_env%ri_rs%grid_cache(iatom)%npts = &
            SIZE(bs_env%ri_rs%grid_cache(iatom)%raw_points, 2)
      END DO
      IF (unit_nr > 0) FLUSH (unit_nr)
 
      bs_env%ri_rs%Z_lP_exists = .false.
      CALL gw_3c_ctx_release(integral_context)
      CALL timestop(handle)
 
   END SUBROUTINE optimize_ri_rs_grid
 
! **************************************************************************************************
!> \brief Evaluate the regularized three-centre-integral fitting error.
!> \param coordinates Flattened atom-relative grid coordinates.
!> \param bs_env GW environment containing RI-RS configuration and grid data.
!> \param clusters Rank-local clusters and their exact three-centre integrals.
!> \param grid_offsets Starting coordinate offset for each atom.
!> \param value Mean normalized squared three-centre-integral error.
!> \param gradient Derivative of value with respect to coordinates.
!> \param maximum_absolute_error Largest absolute three-centre-integral error.
!> \param valid Whether every local-cluster evaluation succeeded.
! **************************************************************************************************
   SUBROUTINE grid_objective(coordinates, bs_env, clusters, grid_offsets, value, gradient, &
                             maximum_absolute_error, valid)
      REAL(kind=dp), DIMENSION(:), INTENT(IN)            :: coordinates
      TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
      TYPE(local_cluster_type), DIMENSION(:), INTENT(IN) :: clusters
      INTEGER, DIMENSION(:), INTENT(IN)                  :: grid_offsets
      REAL(kind=dp), INTENT(OUT)                         :: value
      REAL(kind=dp), DIMENSION(:), INTENT(OUT)           :: gradient
      REAL(kind=dp), INTENT(OUT)                         :: maximum_absolute_error
      LOGICAL, INTENT(OUT)                               :: valid
 
      INTEGER                                            :: all_valid, icluster, successful_clusters
      REAL(kind=dp)                                      :: cluster_maximum_absolute_error, &
                                                            cluster_value
      REAL(kind=dp), ALLOCATABLE                         :: cluster_gradient(:, :)
      TYPE(cell_type), POINTER                           :: cell
      TYPE(mp_para_env_type), POINTER                    :: para_env
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
 
      cell => bs_env%ri_rs%cell
      particle_set => bs_env%ri_rs%particle_set
      para_env => bs_env%para_env
      CALL unpack_atom_grids(coordinates, bs_env%ri_rs%grid_cache, grid_offsets)
      value = 0.0_dp
      gradient = 0.0_dp
      successful_clusters = 0
      maximum_absolute_error = 0.0_dp
      valid = .true.
      DO icluster = 1, SIZE(clusters)
         CALL evaluate_local_cluster(clusters(icluster), bs_env%ri_rs%grid_cache, &
                                     bs_env, particle_set, cell, &
                                     cluster_value, cluster_gradient, &
                                     cluster_maximum_absolute_error, valid)
         IF (.NOT. valid) THEN
            IF (ALLOCATED(cluster_gradient)) DEALLOCATE (cluster_gradient)
            EXIT
         END IF
         value = value + cluster_value
         CALL accumulate_atom_gradient(clusters(icluster)%atoms, bs_env%ri_rs%grid_cache, &
                                       grid_offsets, cluster_gradient, gradient)
         successful_clusters = successful_clusters + 1
         maximum_absolute_error = &
            max(maximum_absolute_error, cluster_maximum_absolute_error)
         DEALLOCATE (cluster_gradient)
      END DO
      all_valid = merge(1, 0, valid)
      CALL para_env%sum(value)
      CALL para_env%sum(gradient)
      CALL para_env%sum(successful_clusters)
      CALL para_env%sum(all_valid)
      CALL para_env%max(maximum_absolute_error)
      valid = successful_clusters == SIZE(particle_set) .AND. all_valid == para_env%num_pe
      IF (.NOT. valid) THEN
         value = huge(value)
         gradient = 0.0_dp
         RETURN
      END IF
      value = value/real(SIZE(particle_set), dp)
      gradient = gradient/real(SIZE(particle_set), dp)
   END SUBROUTINE grid_objective
 
! **************************************************************************************************
!> \brief Minimize the RI-RS grid objective with CP2K's bound-constrained L-BFGS implementation.
!> \param x Coordinates on entry and best coordinates found on return.
!> \param lower Lower bound for each coordinate.
!> \param upper Upper bound for each coordinate.
!> \param accuracy Projected-gradient convergence threshold.
!> \param max_evaluations Maximum number of objective evaluations.
!> \param bs_env GW environment containing RI-RS configuration and grid data.
!> \param clusters Rank-local clusters and their exact three-centre integrals.
!> \param grid_offsets Starting coordinate offset for each atom.
! **************************************************************************************************
   SUBROUTINE optimize_grid_coordinates(x, lower, upper, accuracy, max_evaluations, &
                                        bs_env, clusters, grid_offsets)
      REAL(kind=dp), DIMENSION(:), INTENT(INOUT)         :: x
      REAL(kind=dp), DIMENSION(:), INTENT(IN)            :: lower, upper
      REAL(kind=dp), INTENT(IN)                          :: accuracy
      INTEGER, INTENT(IN)                                :: max_evaluations
      TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
      TYPE(local_cluster_type), DIMENSION(:), INTENT(IN) :: clusters
      INTEGER, DIMENSION(:), INTENT(IN)                  :: grid_offsets
 
      INTEGER, PARAMETER                                 :: memory = 7
 
      CHARACTER(LEN=60)                                  :: csave, task
      INTEGER                                            :: evaluations
      REAL(kind=dp), DIMENSION(29)                       :: dsave
      REAL(kind=dp), DIMENSION(2*memory*SIZE(x)+5*SIZE(x&
         )+11*memory**2+8*memory)                        :: wa
      REAL(kind=dp)                                      :: best_f, f, maximum_absolute_error
      LOGICAL, DIMENSION(4)                              :: lsave
      LOGICAL                                            :: evaluation_ok, have_best
      INTEGER, DIMENSION(SIZE(x))                        :: bound_type
      INTEGER, DIMENSION(44)                             :: isave
      INTEGER, DIMENSION(3*SIZE(x))                      :: iwa
      REAL(kind=dp), DIMENSION(SIZE(x))                  :: best_x, g
 
      cpassert(SIZE(x) > 0)
      cpassert(all(shape(lower) == shape(x)) .AND. all(shape(upper) == shape(x)))
      cpassert(all(lower <= upper))
      cpassert(accuracy > 0.0_dp .AND. max_evaluations > 0)
 
      bound_type = 2
      task = 'START'
      csave = ''
      f = huge(f)
      g = 0.0_dp
      wa = 0.0_dp
      iwa = 0
      lsave = .false.
      isave = 0
      dsave = 0.0_dp
      evaluations = 0
      have_best = .false.
      best_f = huge(best_f)
      best_x = x
 
      DO
         CALL setulb(SIZE(x), memory, x, lower, upper, bound_type, f, g, &
                     0.0_dp, accuracy, wa, iwa, task, -1, csave, lsave, isave, dsave, -1.0_dp)
         IF (task(1:2) == 'FG') THEN
            IF (evaluations >= max_evaluations) EXIT
            CALL grid_objective(x, bs_env, clusters, grid_offsets, f, g, &
                                maximum_absolute_error, evaluation_ok)
            evaluations = evaluations + 1
            IF (.NOT. evaluation_ok) EXIT
            IF (.NOT. have_best .OR. f < best_f) THEN
               have_best = .true.
               best_f = f
               best_x = x
            END IF
         ELSE IF (task(1:5) == 'NEW_X') THEN
            cycle
         ELSE
            EXIT
         END IF
      END DO
      IF (have_best) x = best_x
   END SUBROUTINE optimize_grid_coordinates
 
! **************************************************************************************************
!> \brief Build local clusters and precompute exact (mu nu|P) values on their owning rank.
!> \param context Three-centre-integral evaluation context.
!> \param particle_set Atomic positions and kinds.
!> \param cell Simulation cell.
!> \param ao_size Number of orbital basis functions on each atom.
!> \param ri_size Number of auxiliary basis functions on each atom.
!> \param cutoff Radius of each atom-centred local cluster.
!> \param para_env MPI environment used to distribute cluster ownership.
!> \param clusters Rank-local clusters and their exact three-centre integrals.
! **************************************************************************************************
   SUBROUTINE build_local_clusters(context, particle_set, cell, ao_size, ri_size, cutoff, &
                                   para_env, clusters)
      TYPE(gw_3c_ctx_type), INTENT(IN)                   :: context
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      TYPE(cell_type), POINTER                           :: cell
      INTEGER, DIMENSION(:), INTENT(IN)                  :: ao_size, ri_size
      REAL(kind=dp), INTENT(IN)                          :: cutoff
      TYPE(mp_para_env_type), POINTER                    :: para_env
      TYPE(local_cluster_type), ALLOCATABLE, INTENT(OUT) :: clusters(:)
 
      INTEGER                                            :: center, ia, iatom, ja, jatom, ka, katom, &
                                                            n_ao, n_cluster, n_ri
      INTEGER, ALLOCATABLE                               :: ao_offset(:), atom_buffer(:), &
                                                            ri_offset(:)
      LOGICAL                                            :: screened
      TYPE(gw_3c_ws_type)                                :: workspace
 
      n_cluster = count([(mod(center - 1, para_env%num_pe) == para_env%mepos, &
                          center=1, SIZE(particle_set))])
      ALLOCATE (clusters(n_cluster), atom_buffer(SIZE(particle_set)), &
                ao_offset(SIZE(particle_set)), ri_offset(SIZE(particle_set)))
      CALL gw_3c_ws_create(workspace, context)
      n_cluster = 0
      DO center = 1, SIZE(particle_set)
         IF (mod(center - 1, para_env%num_pe) /= para_env%mepos) cycle
         n_cluster = n_cluster + 1
         ia = 0
         DO iatom = 1, SIZE(particle_set)
            IF (sum(pbc(particle_set(iatom)%r - particle_set(center)%r, cell)**2) <= cutoff**2) THEN
               ia = ia + 1
               atom_buffer(ia) = iatom
            END IF
         END DO
         ALLOCATE (clusters(n_cluster)%atoms(ia), source=atom_buffer(1:ia))
         ao_offset = 0
         ri_offset = 0
         n_ao = 0
         n_ri = 0
         DO ia = 1, SIZE(clusters(n_cluster)%atoms)
            iatom = clusters(n_cluster)%atoms(ia)
            ao_offset(iatom) = n_ao
            ri_offset(iatom) = n_ri
            n_ao = n_ao + ao_size(iatom)
            n_ri = n_ri + ri_size(iatom)
         END DO
         ALLOCATE (clusters(n_cluster)%three_center(n_ao, n_ao, n_ri), source=0.0_dp)
         DO ia = 1, SIZE(clusters(n_cluster)%atoms)
            iatom = clusters(n_cluster)%atoms(ia)
            DO ka = 1, SIZE(clusters(n_cluster)%atoms)
               katom = clusters(n_cluster)%atoms(ka)
               DO ja = 1, SIZE(clusters(n_cluster)%atoms)
                  jatom = clusters(n_cluster)%atoms(ja)
                  CALL build_3c_integral_block_ctx( &
                     clusters(n_cluster)%three_center, context, workspace, &
                     atom_j=jatom, atom_k=katom, atom_i=iatom, &
                     j_offset=ao_offset(jatom), k_offset=ao_offset(katom), &
                     i_offset=ri_offset(iatom), screened=screened)
               END DO
            END DO
         END DO
      END DO
      CALL gw_3c_ws_release(workspace)
   END SUBROUTINE build_local_clusters
 
! **************************************************************************************************
!> \brief Evaluate one complete cluster and return derivatives for its physical grid points.
!> \param cluster Local atoms and exact three-centre integrals.
!> \param grids Atom-centred RI-RS grids.
!> \param bs_env GW environment containing orbital basis sets and RI-RS parameters.
!> \param particle_set Atomic positions and kinds.
!> \param cell Simulation cell.
!> \param value Normalized squared three-centre-integral error.
!> \param gradient Derivative with respect to the cluster's physical grid coordinates.
!> \param maximum_absolute_error Largest absolute three-centre-integral error.
!> \param successful Whether the regularized fitting equations were solved.
! **************************************************************************************************
   SUBROUTINE evaluate_local_cluster(cluster, grids, bs_env, particle_set, cell, &
                                     value, gradient, maximum_absolute_error, successful)
      TYPE(local_cluster_type), INTENT(IN)               :: cluster
      TYPE(rirs_grid_type), DIMENSION(:), INTENT(IN)     :: grids
      TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      TYPE(cell_type), POINTER                           :: cell
      REAL(kind=dp), INTENT(OUT)                         :: value
      REAL(kind=dp), ALLOCATABLE, INTENT(OUT)            :: gradient(:, :)
      REAL(kind=dp), INTENT(OUT)                         :: maximum_absolute_error
      LOGICAL, INTENT(OUT)                               :: successful
 
      INTEGER                                            :: ao_offset, ia, iatom, ikind, n_ao, &
                                                            n_atom_ao, n_grid, point_offset
      REAL(kind=dp), ALLOCATABLE                         :: dphi(:, :, :), phi(:, :), points(:, :)
      TYPE(gto_basis_set_type), POINTER                  :: basis
 
      n_ao = SIZE(cluster%three_center, 1)
      n_grid = 0
      DO ia = 1, SIZE(cluster%atoms)
         iatom = cluster%atoms(ia)
         n_grid = n_grid + SIZE(grids(iatom)%raw_points, 2)
      END DO
      ALLOCATE (points(3, n_grid), phi(n_grid, n_ao), &
                dphi(3, n_grid, n_ao), gradient(3, n_grid))
      phi = 0.0_dp
      dphi = 0.0_dp
      point_offset = 0
      DO ia = 1, SIZE(cluster%atoms)
         iatom = cluster%atoms(ia)
         points(:, point_offset + 1:point_offset + SIZE(grids(iatom)%raw_points, 2)) = &
            spread(particle_set(iatom)%r, 2, SIZE(grids(iatom)%raw_points, 2)) + &
            grids(iatom)%raw_points
         point_offset = point_offset + SIZE(grids(iatom)%raw_points, 2)
      END DO
      ao_offset = 0
      DO ia = 1, SIZE(cluster%atoms)
         iatom = cluster%atoms(ia)
         ikind = particle_set(iatom)%atomic_kind%kind_number
         basis => bs_env%basis_set_AO(ikind)%gto_basis_set
         n_atom_ao = basis%nsgf
         CALL evaluate_ao_basis_on_points( &
            phi(:, ao_offset + 1:ao_offset + n_atom_ao), points, basis, &
            particle_set(iatom)%r, cell, &
            dphi=dphi(:, :, ao_offset + 1:ao_offset + n_atom_ao))
         ao_offset = ao_offset + n_atom_ao
      END DO
      CALL evaluate_rirs_grid_cluster(phi, dphi, cluster%three_center, bs_env%ri_rs%tikhonov, &
                                      value, gradient, maximum_absolute_error, successful)
   END SUBROUTINE evaluate_local_cluster
 
! **************************************************************************************************
!> \brief Evaluate the regularized local-cluster three-centre fit and its coordinate gradient.
!>
!> The solve applies the same column Jacobi scaling and Tikhonov parameter as the production
!> Z_lP construction. The objective contains only the normalized three-centre-integral residual.
!> \param phi AO values, indexed (l,mu).
!> \param dphi Cartesian derivatives of AO values, indexed (alpha,l,mu).
!> \param three_center Exact three-centre integrals, indexed (mu,nu,P).
!> \param tikhonov Tikhonov parameter used by the production RI-RS solve.
!> \param value Normalized squared residual for this cluster.
!> \param gradient Analytic derivative d value/d r(alpha,l).
!> \param maximum_absolute_error Largest absolute error in an unweighted three-centre integral.
!> \param successful False if the regularized normal equations cannot be solved.
! **************************************************************************************************
   SUBROUTINE evaluate_rirs_grid_cluster(phi, dphi, three_center, tikhonov, value, gradient, &
                                         maximum_absolute_error, successful)
      REAL(kind=dp), DIMENSION(:, :), INTENT(IN)         :: phi
      REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN)      :: dphi, three_center
      REAL(kind=dp), INTENT(IN)                          :: tikhonov
      REAL(kind=dp), INTENT(OUT)                         :: value
      REAL(kind=dp), DIMENSION(:, :), INTENT(OUT)        :: gradient
      REAL(kind=dp), INTENT(OUT)                         :: maximum_absolute_error
      LOGICAL, INTENT(OUT)                               :: successful
 
      CHARACTER(len=*), PARAMETER :: routinen = 'evaluate_rirs_grid_cluster'
      REAL(kind=dp), PARAMETER                           :: jacobi_floor = 1.0e-16_dp
 
      INTEGER                                            :: alpha, handle, info, ipair, l, mu, n_ao, &
                                                            n_grid, n_pair, n_ri, nu, p
      INTEGER, ALLOCATABLE                               :: pair_mu(:), pair_nu(:)
      REAL(kind=dp)                                      :: absolute_error, column_norm2, d_a, &
                                                            denom, factor, projection, scale
      REAL(kind=dp), ALLOCATABLE :: a_matrix(:, :), a_scaled(:, :), b_matrix(:, :), &
         derivative_a(:, :), derivative_scaled_a(:), gram(:, :), residual(:, :), rhs(:, :), &
         scale_columns(:), y_matrix(:, :), yz_matrix(:, :), z_matrix(:, :)
 
      CALL timeset(routinen, handle)
 
      n_grid = SIZE(phi, 1)
      n_ao = SIZE(phi, 2)
      n_ri = SIZE(three_center, 3)
      n_pair = n_ao*(n_ao + 1)/2
 
      cpassert(SIZE(dphi, 1) == 3)
      cpassert(SIZE(dphi, 2) == n_grid)
      cpassert(SIZE(dphi, 3) == n_ao)
      cpassert(SIZE(three_center, 1) == n_ao)
      cpassert(SIZE(three_center, 2) == n_ao)
      cpassert(SIZE(gradient, 1) == 3)
      cpassert(SIZE(gradient, 2) == n_grid)
      cpassert(tikhonov >= 0.0_dp)
 
      value = huge(value)
      gradient = 0.0_dp
      maximum_absolute_error = huge(maximum_absolute_error)
      successful = .false.
      IF (n_grid < 1 .OR. n_ri < 1) THEN
         CALL timestop(handle)
         RETURN
      END IF
 
      ALLOCATE (a_matrix(n_pair, n_grid), a_scaled(n_pair, n_grid), &
                b_matrix(n_pair, n_ri), pair_mu(n_pair), pair_nu(n_pair), &
                scale_columns(n_grid))
      ipair = 0
      DO nu = 1, n_ao
         DO mu = 1, nu
            ipair = ipair + 1
            pair_mu(ipair) = mu
            pair_nu(ipair) = nu
            factor = sqrt(real(2 - merge(1, 0, mu == nu), dp))
            DO l = 1, n_grid
               a_matrix(ipair, l) = factor*phi(l, mu)*phi(l, nu)
            END DO
            DO p = 1, n_ri
               b_matrix(ipair, p) = factor*three_center(mu, nu, p)
            END DO
         END DO
      END DO
 
      denom = sum(b_matrix*b_matrix)
      IF (denom <= tiny(1.0_dp)) THEN
         DEALLOCATE (a_matrix, a_scaled, b_matrix, pair_mu, pair_nu, scale_columns)
         CALL timestop(handle)
         RETURN
      END IF
 
      ! Production Z_lP uses D'=d(A^T A)d+lambda I with d_l=1/sqrt((A^T A)_ll).
      DO l = 1, n_grid
         column_norm2 = sum(a_matrix(:, l)*a_matrix(:, l))
         scale_columns(l) = 1.0_dp/sqrt(max(column_norm2, jacobi_floor))
         a_scaled(:, l) = scale_columns(l)*a_matrix(:, l)
      END DO
      ALLOCATE (gram(n_grid, n_grid), rhs(n_grid, n_ri), z_matrix(n_grid, n_ri))
      CALL dgemm('T', 'N', n_grid, n_grid, n_pair, 1.0_dp, a_scaled, n_pair, &
                 a_scaled, n_pair, 0.0_dp, gram, n_grid)
      DO l = 1, n_grid
         gram(l, l) = gram(l, l) + tikhonov
      END DO
      CALL dgemm('T', 'N', n_grid, n_ri, n_pair, 1.0_dp, a_scaled, n_pair, &
                 b_matrix, n_pair, 0.0_dp, rhs, n_grid)
      CALL dpotrf('L', n_grid, gram, n_grid, info)
      IF (info /= 0) THEN
         DEALLOCATE (a_matrix, a_scaled, b_matrix, gram, pair_mu, pair_nu, rhs, &
                     scale_columns, z_matrix)
         CALL timestop(handle)
         RETURN
      END IF
      z_matrix(:, :) = rhs
      CALL dpotrs('L', n_grid, n_ri, gram, n_grid, z_matrix, n_grid, info)
      DEALLOCATE (rhs)
      IF (info /= 0) THEN
         DEALLOCATE (a_matrix, a_scaled, b_matrix, gram, pair_mu, pair_nu, &
                     scale_columns, z_matrix)
         CALL timestop(handle)
         RETURN
      END IF
 
      ALLOCATE (residual(n_pair, n_ri))
      residual(:, :) = b_matrix
      CALL dgemm('N', 'N', n_pair, n_ri, n_grid, -1.0_dp, a_scaled, n_pair, &
                 z_matrix, n_grid, 1.0_dp, residual, n_pair)
      value = sum(residual*residual)/denom
 
      ! For lambda > 0, the residual-only derivative includes the response of the regularized
      ! coefficients. With Y=(A^T A+lambda I)^(-1) Z:
      ! dE/dA = [-2 R Z^T - 2 lambda (R Y^T - A Y Z^T)] / ||B||^2.
      ALLOCATE (y_matrix(n_grid, n_ri), derivative_a(n_pair, n_grid))
      y_matrix(:, :) = z_matrix
      CALL dpotrs('L', n_grid, n_ri, gram, n_grid, y_matrix, n_grid, info)
      IF (info /= 0) THEN
         DEALLOCATE (a_matrix, a_scaled, b_matrix, derivative_a, gram, pair_mu, pair_nu, &
                     residual, scale_columns, y_matrix, z_matrix)
         CALL timestop(handle)
         RETURN
      END IF
      CALL dgemm('N', 'T', n_pair, n_grid, n_ri, -2.0_dp/denom, residual, n_pair, &
                 z_matrix, n_grid, 0.0_dp, derivative_a, n_pair)
      IF (tikhonov > 0.0_dp) THEN
         CALL dgemm('N', 'T', n_pair, n_grid, n_ri, -2.0_dp*tikhonov/denom, residual, n_pair, &
                    y_matrix, n_grid, 1.0_dp, derivative_a, n_pair)
         ALLOCATE (yz_matrix(n_grid, n_grid))
         CALL dgemm('N', 'T', n_grid, n_grid, n_ri, 1.0_dp, y_matrix, n_grid, &
                    z_matrix, n_grid, 0.0_dp, yz_matrix, n_grid)
         CALL dgemm('N', 'N', n_pair, n_grid, n_grid, 2.0_dp*tikhonov/denom, a_scaled, n_pair, &
                    yz_matrix, n_grid, 1.0_dp, derivative_a, n_pair)
         DEALLOCATE (yz_matrix)
      END IF
 
      ALLOCATE (derivative_scaled_a(n_pair))
      DO l = 1, n_grid
         scale = scale_columns(l)
         column_norm2 = sum(a_matrix(:, l)*a_matrix(:, l))
         DO alpha = 1, 3
            DO ipair = 1, n_pair
               mu = pair_mu(ipair)
               nu = pair_nu(ipair)
               factor = sqrt(real(2 - merge(1, 0, mu == nu), dp))
               d_a = factor*(dphi(alpha, l, mu)*phi(l, nu) + &
                             phi(l, mu)*dphi(alpha, l, nu))
               derivative_scaled_a(ipair) = scale*d_a
            END DO
            IF (column_norm2 > jacobi_floor) THEN
               projection = dot_product(a_matrix(:, l), derivative_scaled_a)/scale
               derivative_scaled_a(:) = derivative_scaled_a - &
                                        scale**3*a_matrix(:, l)*projection
            END IF
            gradient(alpha, l) = dot_product(derivative_a(:, l), derivative_scaled_a)
         END DO
      END DO
 
      maximum_absolute_error = 0.0_dp
      DO p = 1, n_ri
         DO ipair = 1, n_pair
            factor = sqrt(real(2 - merge(1, 0, pair_mu(ipair) == pair_nu(ipair)), dp))
            absolute_error = abs(residual(ipair, p))/factor
            maximum_absolute_error = max(maximum_absolute_error, absolute_error)
         END DO
      END DO
 
      successful = .true.
      DEALLOCATE (a_matrix, a_scaled, b_matrix, derivative_a, derivative_scaled_a, gram, &
                  pair_mu, pair_nu, residual, scale_columns, y_matrix, z_matrix)
      CALL timestop(handle)
   END SUBROUTINE evaluate_rirs_grid_cluster
 
! **************************************************************************************************
!> \brief Sum each physical occurrence into its atom-relative coordinate.
!> \param atoms Atom indices in the local cluster.
!> \param grids Atom-centred RI-RS grids.
!> \param grid_offsets Starting coordinate offset for each atom.
!> \param physical_gradient Gradient for the cluster's physical grid points.
!> \param gradient Global flattened atom-relative gradient to update.
! **************************************************************************************************
   SUBROUTINE accumulate_atom_gradient(atoms, grids, grid_offsets, physical_gradient, gradient)
      INTEGER, DIMENSION(:), INTENT(IN)                  :: atoms
      TYPE(rirs_grid_type), DIMENSION(:), INTENT(IN)     :: grids
      INTEGER, DIMENSION(:), INTENT(IN)                  :: grid_offsets
      REAL(kind=dp), DIMENSION(:, :), INTENT(IN)         :: physical_gradient
      REAL(kind=dp), DIMENSION(:), INTENT(INOUT)         :: gradient
 
      INTEGER                                            :: alpha, ia, iatom, l, point_offset
 
      point_offset = 0
      DO ia = 1, SIZE(atoms)
         iatom = atoms(ia)
         DO l = 1, SIZE(grids(iatom)%raw_points, 2)
            DO alpha = 1, 3
               gradient(grid_offsets(iatom) + 3*(l - 1) + alpha) = &
                  gradient(grid_offsets(iatom) + 3*(l - 1) + alpha) + &
                  physical_gradient(alpha, point_offset + l)
            END DO
         END DO
         point_offset = point_offset + SIZE(grids(iatom)%raw_points, 2)
      END DO
   END SUBROUTINE accumulate_atom_gradient
 
! **************************************************************************************************
!> \brief Flatten atom-relative Cartesian coordinates into the optimizer vector.
!> \param grids Atom-centred RI-RS grids.
!> \param grid_offsets Starting coordinate offset for each atom.
!> \param coordinates Flattened atom-relative grid coordinates.
! **************************************************************************************************
   SUBROUTINE pack_atom_grids(grids, grid_offsets, coordinates)
      TYPE(rirs_grid_type), DIMENSION(:), INTENT(IN)     :: grids
      INTEGER, DIMENSION(:), INTENT(IN)                  :: grid_offsets
      REAL(kind=dp), DIMENSION(:), INTENT(OUT)           :: coordinates
 
      INTEGER                                            :: iatom, n_coordinate
 
      cpassert(SIZE(grid_offsets) == SIZE(grids))
      DO iatom = 1, SIZE(grids)
         n_coordinate = SIZE(grids(iatom)%raw_points)
         coordinates(grid_offsets(iatom) + 1:grid_offsets(iatom) + n_coordinate) = &
            reshape(grids(iatom)%raw_points, [n_coordinate])
      END DO
   END SUBROUTINE pack_atom_grids
 
! **************************************************************************************************
!> \brief Restore atom-relative Cartesian coordinates from the optimizer vector.
!> \param coordinates Flattened atom-relative grid coordinates.
!> \param grids Atom-centred RI-RS grids to update.
!> \param grid_offsets Starting coordinate offset for each atom.
! **************************************************************************************************
   SUBROUTINE unpack_atom_grids(coordinates, grids, grid_offsets)
      REAL(kind=dp), DIMENSION(:), INTENT(IN)            :: coordinates
      TYPE(rirs_grid_type), DIMENSION(:), INTENT(INOUT)  :: grids
      INTEGER, DIMENSION(:), INTENT(IN)                  :: grid_offsets
 
      INTEGER                                            :: iatom, n_coordinate
 
      cpassert(SIZE(grid_offsets) == SIZE(grids))
      DO iatom = 1, SIZE(grids)
         n_coordinate = SIZE(grids(iatom)%raw_points)
         grids(iatom)%raw_points(:, :) = reshape( &
                                         coordinates(grid_offsets(iatom) + 1:grid_offsets(iatom) + n_coordinate), &
                                         shape(grids(iatom)%raw_points))
      END DO
   END SUBROUTINE unpack_atom_grids
 
END MODULE gw_optimize_ri_rs_grid