test_isdf.F90

Go to the documentation of this file.

!! Copyright (C) 2024 - 2025 A 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, USA.
 
#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 centroids_oct_m
  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 parser_oct_m
  use profiling_oct_m
  use scf_oct_m
  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_options_oct_m
  use isdf_serial_oct_m
  use isdf_oct_m
  use isdf_utils_oct_m
 
  implicit none
  private
  public :: test_isdf
 
contains
 
  subroutine test_isdf(namespace, serial)
    type(namespace_t),  intent(in) :: namespace
    logical, optional,  intent(in) :: serial
 
    type(electrons_t),     pointer :: sys
    real(real64),        parameter :: LMM_R_SINGLE_ATOM = 100.0_real64
 
    ! Centroids
    type(centroids_t)              :: centroids
    integer(int64),    allocatable :: indices(:)
    integer                        :: n_int
 
    ! ISDF
    integer                        :: n_ks_states
    type(isdf_options_t)           :: isdf_options
    logical                        :: output_cubes, run_serial_isdf
    real(real64),      allocatable :: phi_mu(:, :), P_r_mu(:, :)
    real(real64),      allocatable :: isdf_vectors(:, :), gram_matrix(:, :)
 
    push_sub(test_isdf)
 
    run_serial_isdf = optional_default(serial, .false.)
 
    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 parse_variable(namespace, 'ACESize', sys%st%nst, n_ks_states)
 
    call isdf_options%init(namespace, sys%space, sys%st, sys%gr, n_ks_states)
 
    call test_isdf_get_interpolation_points(namespace, isdf_options, sys, centroids)
 
    ! Visualisation exact and approx product functions
    output_cubes = .false.
 
    if (run_serial_isdf) then
      if (sys%st%parallel_in_states .or. sys%gr%parallel_in_domains) then
        write(message(1), '(A)') 'Reference serial ISDF test requested, but state or domain parallel'
        write(message(2), '(A)') ' specified in the input.'
        call messages_fatal(2, namespace=namespace)
      endif
 
      indices = centroids%global_mesh_indices()
      n_int = size(indices)
 
      call isdf_serial_interpolation_vectors(isdf_options, namespace, sys%gr, sys%st, indices, &
        phi_mu, p_r_mu, isdf_vectors)
      safe_deallocate_a(phi_mu)
      safe_deallocate_a(p_r_mu)
 
      safe_allocate(gram_matrix(1:n_int, 1:n_int))
      call isdf_gram_matrix(sys%gr, isdf_vectors, gram_matrix)
      call output_matrix(namespace, "isdf_gram.dat", gram_matrix)
 
      ! This deallocates isdf_vectors to save memory
      call quantify_error_and_visualise(isdf_options, namespace, sys%st, sys%space, sys%gr, &
        sys%ions, indices, isdf_vectors, output_cubes)
 
      deallocate(indices)
 
    else
      call messages_write('Info: Testing ISDF Parallel Implementation')
      call messages_new_line()
      call messages_info(namespace=namespace)
 
      n_int = centroids%npoints_global()
      if (mpi_world%is_root()) write(*, *) 'ISDF: Total number of centroids: ', n_int
      write(*, *) "ISDF: Number of centroids local to domain process ", sys%gr%mpi_grp%rank, " is ", centroids%npoints()
 
      call isdf_interpolation_vectors(isdf_options, namespace, sys%gr, sys%st, centroids, phi_mu, p_r_mu, isdf_vectors)
      safe_deallocate_a(phi_mu)
      safe_deallocate_a(p_r_mu)
 
      ! Gram matrix of the interpolation vectors offers a quantity reduced over states and domains
      ! so good for testing parallelisation on
      safe_allocate(gram_matrix(1:n_int, 1:n_int))
      call isdf_gram_matrix(sys%gr, isdf_vectors, gram_matrix)
      call output_matrix(namespace, "isdf_gram.dat", gram_matrix)
    endif
 
    ! Tear down
    safe_deallocate_a(gram_matrix)
    safe_deallocate_a(isdf_vectors)
    call states_elec_deallocate_wfns(sys%st)
    safe_deallocate_p(sys)
 
    pop_sub(test_isdf)
 
  end subroutine test_isdf
 
 
  subroutine test_isdf_get_interpolation_points(namespace, isdf_options, sys, centroids)
    type(namespace_t),               intent(in   ) :: namespace
    type(isdf_options_t),            intent(inout) :: isdf_options
    type(electrons_t),  pointer,     intent(in   ) :: sys
 
    type(centroids_t),               intent(out) :: centroids
 
    character(len=1),  allocatable :: dummy_species(:)
 
    push_sub(test_isdf_get_interpolation_points)
 
    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",  isdf_options%interpolation_coords, dummy_species)
    endif
 
    call weighted_kmeans(sys%space, sys%gr,  sys%st%rho(1:sys%gr%np, 1), isdf_options%interpolation_coords)
 
    if (debug%info) then
      dummy_species = 'B'
      call write_standard_xyz_file(namespace, static_dir // "/final_centroids",  isdf_options%interpolation_coords, dummy_species)
    endif
 
    ! Convert positions to indices
    call centroids%init(sys%gr, isdf_options%interpolation_coords, check_duplicates=.true.)
 
    ! Output the centroid indices so they can be asserted on
    call centroids%output_all_indices(namespace, sys%grp, "centroid_indices.dat")
 
    pop_sub(test_isdf_get_interpolation_points)
 
  end subroutine test_isdf_get_interpolation_points
 
end module test_isdf_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: