isdf_options.F90

Go to the documentation of this file.

!! Copyright (C) 2025. Alexander Buccheri
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, U
 
#include "global.h"
 
module isdf_options_oct_m
  use, intrinsic :: iso_fortran_env
 
  use centroids_oct_m,         only: centroids_t
  use debug_oct_m
  use global_oct_m
  use io_oct_m
  use kmeans_clustering_oct_m, only: weighted_kmeans, sample_initial_centroids
  use mesh_oct_m
  use messages_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use space_oct_m
  use states_elec_oct_m
  implicit none
  private
 
  type, public :: isdf_options_t
    logical                     :: check_n_interp
    integer                     :: n_interp
    integer                     :: n_ks_states = 0
    !                                                     This should be consistent with the subspace of ACE
    integer(int64)              :: seed
    logical                     :: use_serial
    type(centroids_t)           :: centroids
    real(real64),   allocatable :: interpolation_coords(:, :)
 
  contains
    procedure :: init => isdf_options_init
    procedure :: write_info => isdf_options_write_info
    procedure :: get_interpolation_points => isdf_get_interpolation_points
    procedure :: end => isdf_options_end
    final     :: isdf_options_finalise
  end type isdf_options_t
 
  character(len=0), parameter, private :: NULL_FILE = ""
 
contains
 
  subroutine isdf_options_init(this, namespace, space, gr, n_ks_states, rng_seed)
    class(isdf_options_t),     intent(inout) :: this
    type(namespace_t),         intent(in )   :: namespace
    class(space_t),            intent(in )   :: space
    class(mesh_t),             intent(in )   :: gr
    integer,                   intent(in )   :: n_ks_states
    ! If ISDF is used with ACE, the number of states should be greater or equal to the subspace in which ACE
    ! is exact. Else the density fitting will not capture the features of these states.
    integer(int64),            intent(inout), optional :: rng_seed
 
    integer,        parameter   :: NULL_N_INTERP = 0
    integer(int64), parameter   :: default_rng_seed = 101_int64
 
    push_sub(isdf_options_init)
 
    this%seed =  optional_default(rng_seed, default_rng_seed)
    this%n_ks_states = n_ks_states
 
    !%Variable ISDFNpoints
    !%Type integer
    !%Default 0
    !%Section ISDF
    !%Description
    !% Specify the number of interpolation points, and therefore the number density fitting vectors.
    !% The number of interpolation points should be considerably smaller than the number
    !% of grid points. ISDFNpoints MUST be set. A sensible guess is 10 times the number of occupied states.
    !%End
 
    !%Variable CheckISDFNpoints
    !%Type logical
    !%Default no
    !%Section ISDF
    !%Description
    !% We can get a very good indication of the optimal number of interpolation points to use in ISDF
    !% by choosing an excessive number,and computing the rank of the fitting coefficient Gram matrix, $CC^T$.
    !% If the rank is equal to ISDFNpoints, this suggests one can choose more, or has the optimal number.
    !% If the rank is less than ISDFNpoints, then it defines the upper bound beyond which the exact exchange energy
    !% will not increase in precision. The same argument is not true of the two-electron integrals, which can be
    !% converged further at the expense of an increase in the condition number of $CC^T$.
    !%End
 
    !%Variable ISDFRunSerial
    !%Type logical
    !%Default no
    !%Section ISDF
    !%Description
    !% Use the reference serial implementation of ISDF.
    !%
    !%End
    call parse_variable(namespace, 'ISDFRunSerial', .false., this%use_serial)
    call parse_variable(namespace, 'CheckISDFNpoints', .false., this%check_n_interp)
 
    ! Set number of interpolation points
    call parse_variable(namespace, 'ISDFNpoints', null_n_interp, this%n_interp)
    if (this%n_interp == null_n_interp) then
      message(1) = "ISDFNpoints must be specified"
      call messages_fatal(1, namespace=namespace)
    endif
 
    safe_allocate(this%interpolation_coords(1:space%dim, 1:this%n_interp))
    call sample_initial_centroids(gr, this%interpolation_coords, seed_value=this%seed)
 
    pop_sub(isdf_options_init)
 
  end subroutine isdf_options_init
 
 
  subroutine isdf_options_end(this)
    class(isdf_options_t) :: this
 
    push_sub(isdf_options_end)
 
    safe_deallocate_a(this%interpolation_coords)
    call this%centroids%end()
 
    pop_sub(isdf_options_end)
 
  end subroutine isdf_options_end
 
 
  subroutine isdf_options_finalise(this)
    type(isdf_options_t) :: this
 
    push_sub(isdf_options_finalise)
    call this%end()
    pop_sub(isdf_options_finalise)
 
  end subroutine isdf_options_finalise
 
 
  subroutine isdf_options_write_info(this, namespace)
    class(isdf_options_t), intent(in) :: this
    type(namespace_t),    intent(in) :: namespace
 
    push_sub(isdf_options_write_info)
 
    call messages_print_var_value("Interpolation Point Algorithm", "Weighted Kmeans Clustering", &
      namespace=namespace)
 
    call messages_print_var_value("Number of Interpolation Points (ISDFNpoints)", this%n_interp, &
      namespace=namespace)
 
    call messages_print_var_value("Number of states used in the ISDF vector construction", &
      this%n_ks_states, namespace=namespace)
 
    pop_sub(isdf_options_write_info)
 
  end subroutine isdf_options_write_info
 
 
  subroutine isdf_get_interpolation_points(this, namespace, space, mesh, weight)
    class(isdf_options_t),     intent(inout) :: this
    type(namespace_t),         intent(in   ) :: namespace
    class(space_t),            intent(in   ) :: space
    class(mesh_t),             intent(in   ) :: mesh
    real(real64),              intent(in   ) :: weight(:)
 
    logical                      :: check_duplicates
 
    push_sub_with_profile(isdf_get_interpolation_points)
 
    check_duplicates = debug%info .eqv. .true.
 
    write(message(1),'(a, I3)') 'Info: Performing weighted KMeans'
    call messages_info(1, namespace=namespace)
    call weighted_kmeans(space, mesh, weight, this%interpolation_coords)
    call this%centroids%end()
    call this%centroids%init(mesh, this%interpolation_coords, check_duplicates = check_duplicates)
 
    pop_sub_with_profile(isdf_get_interpolation_points)
 
  end subroutine isdf_get_interpolation_points
 
end module isdf_options_oct_m