test_isdf.F90

Go to the documentation of this file.

!! Copyright (C) 2024 A Buccheri.
 
#include "global.h"
 
module test_isdf_oct_m
  use, intrinsic :: iso_fortran_env
 
  use calc_mode_par_oct_m,   only: calc_mode_par, p_strategy_states
  use comm_oct_m
  use debug_oct_m
  use electrons_oct_m,       only: electrons_t
  use face_splitting_oct_m
  use global_oct_m
  use io_function_oct_m
  use io_oct_m,              only: io_open, io_close, io_mkdir
  use lalg_adv_oct_m
  use lcao_oct_m,            only: lcao_run
  use mesh_oct_m
  use messages_oct_m
  use mpi_oct_m,             only: mpi_world, mpi_double_precision
  use namespace_oct_m,       only: namespace_t
  use profiling_oct_m
  use scf_oct_m
  use sort_oct_m,            only: sort
  use states_elec_oct_m
  use types_oct_m
  use unit_system_oct_m,     only: unit_one, unit_angstrom
  use wfs_elec_oct_m
 
  ! Modules under test
  use kmeans_clustering_oct_m
  use isdf_serial_oct_m
 
  implicit none
  private
  public :: test_isdf
 
contains
 
  subroutine map_centroid_positions_to_indices(mesh, positions, indices)
    class(mesh_t),    intent(in ) :: mesh
    real(real64),     intent(in ) :: positions(:, :)
    integer(int64),   intent(out) :: indices(:)
 
    integer(int64) :: ipg
    integer        :: n_centroids, ic
    integer        :: ix(size(positions, 1))
    real(real64)   :: chi(size(positions, 1))
 
    push_sub(map_centroid_positions_to_indices)
 
    if (mesh%parallel_in_domains) then
      message(1) = "map_centroid_positions_to_indices is not domain-parallel"
      call messages_fatal(1)
    endif
 
    n_centroids = size(positions, 2)
    assert(size(indices) == n_centroids)
 
    do ic = 1, n_centroids
      chi = mesh%coord_system%from_cartesian(positions(:, ic))
      ix = nint(chi / mesh%spacing)
      ipg = mesh_global_index_from_coords(mesh, ix)
      indices(ic) = ipg
    enddo
 
    pop_sub(map_centroid_positions_to_indices)
 
  end subroutine map_centroid_positions_to_indices
 
 
  subroutine test_isdf(namespace)
    type(namespace_t),  intent(in) :: namespace
 
    type(electrons_t),     pointer :: sys
    integer                        :: i
    real(real64),        parameter :: LMM_R_SINGLE_ATOM = 100.0_real64
 
    ! Centroids
    integer                        :: n_duplicates
    integer                        :: cnt, unit
    integer(int64)                 :: rng_seed
    integer(int64),    allocatable :: indices(:), unique_indices(:)
    real(real64),      allocatable :: centroid_positions(:, :)
    character(len=1),  allocatable :: dummy_species(:)
 
    ! ISDF
    type(isdf_options_t)           :: isdf_options
    logical                        :: output_cubes
    real(real64),      allocatable :: phi(:, :)
    real(real64),      allocatable :: isdf_vectors(:, :)
 
    push_sub(test_isdf)
 
    call messages_write('Info: Testing ISDF')
    call messages_new_line()
    call messages_new_line()
    call messages_info(namespace=namespace)
 
    ! Set up a system of electrons
    call calc_mode_par%set_parallelization(p_strategy_states, default=.false.)
    ! If generate_epot is false, the code will crash in v_ks_calc
    sys => electrons_t(namespace, generate_epot=.true.)
    call sys%init_parallelization(mpi_world)
 
    ! Initialise a view of the wave functions
    ! Note, this does not change the underlying data layout of the states
    call states_elec_allocate_wfns(sys%st, sys%gr, wfs_type = type_float, packed=.true.)
 
    ! Guess at the states using lcao
    call lcao_run(namespace, sys%space, sys%gr, sys%ions, sys%ext_partners, &
      sys%st, sys%ks, sys%hm, lmm_r = lmm_r_single_atom)
 
    call isdf_options%init(namespace, sys%st)
 
    ! Perform kmeans - Random seed fixed for test reproducibility
    rng_seed = int(101, int64)
    safe_allocate(centroid_positions(1:sys%space%dim, 1:isdf_options%n_interp))
 
    call sample_initial_centroids(sys%gr, centroid_positions, seed_value=rng_seed)
    if (debug%info) then
      allocate(dummy_species(isdf_options%n_interp), source='A')
      call io_mkdir(static_dir, namespace)
      call write_standard_xyz_file(namespace, static_dir // "/initial_centroids",  centroid_positions, dummy_species)
    endif
 
    call weighted_kmeans(sys%space, sys%gr,  sys%st%rho(1:sys%gr%np, 1), centroid_positions)
 
    if (debug%info) then
      dummy_species = 'B'
      call write_standard_xyz_file(namespace, static_dir // "/final_centroids",  centroid_positions, dummy_species)
    endif
 
    ! Convert positions to indices
    safe_allocate(indices(1:isdf_options%n_interp))
    call map_centroid_positions_to_indices(sys%gr, centroid_positions, indices)
 
    ! Output the centroid indices so they can be asserted on
    if (mpi_world%is_root()) then
      open(newunit=unit, file="centroid_indices.dat")
      write(unit, *) 'N centroids', isdf_options%n_interp
      do i = 1, isdf_options%n_interp
        write(unit, *) indices(i)
      enddo
      close(unit)
    endif
 
    ! Confirm no duplicates are present
    call sort(indices)
 
    ! For the sake of testing, remove any duplicates
    allocate(unique_indices(isdf_options%n_interp))
 
    n_duplicates = 0
    cnt = 1
    unique_indices(1) = indices(1)
    do i = 1, isdf_options%n_interp - 1
      if (indices(i) == indices(i+1)) then
        n_duplicates = n_duplicates + 1
        cycle
      endif
      cnt = cnt + 1
      unique_indices(cnt) = indices(i+1)
    enddo
    safe_deallocate_a(indices)
    allocate(indices(cnt), source=unique_indices(1:cnt))
    deallocate(unique_indices)
 
    if (n_duplicates > 0) then
      write(message(1), '(I0)') n_duplicates
      message(1) = trim(message(1)) // " duplicates found after discretisation of centroid positions"
      write(message(2), '(I0)') cnt
      message(2) = " Removed the duplicate/s and resized 'indices' array to " // trim(message(2))
      call messages_warning(2)
      !call messages_fatal(1)
    endif
 
    ! ISDF Testing - Serial code
    if (.not. sys%st%parallel_in_states .and. .not. sys%gr%parallel_in_domains) then
 
      ! Returns isdf_vectors, and phi for reuse
      call serial_interpolative_separable_density_fitting_vectors(namespace, &
        sys%gr, sys%st, indices, phi, isdf_vectors)
 
      ! Error quantification and visualisation
      output_cubes = .false.
 
      ! This deallocates phi and isdf_vectors to save memory
      call quantify_error_and_visualise(namespace, sys%st, sys%space, sys%gr, &
        sys%ions, phi, indices, isdf_vectors, output_cubes)
    endif
 
    call states_elec_deallocate_wfns(sys%st)
    safe_deallocate_p(sys)
    safe_deallocate_a(indices)
 
    pop_sub(test_isdf)
 
  end subroutine test_isdf
 
end module test_isdf_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: