test_kmeans_clustering.F90

Go to the documentation of this file.

!! Copyright (C) 2024 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_kmeans_clustering_oct_m
  use, intrinsic :: iso_fortran_env
  use calc_mode_par_oct_m,   only: calc_mode_par, p_strategy_states
  use debug_oct_m
  use electrons_oct_m,       only: electrons_t
  use global_oct_m
  use io_function_oct_m
  use io_oct_m,              only: io_open, io_close, io_mkdir
  use lcao_oct_m,            only: lcao_run
  use messages_oct_m
  use mpi_oct_m,             only: mpi_world
  use namespace_oct_m,       only: namespace_t
  use profiling_oct_m
  use states_elec_oct_m
  use test_parameters_oct_m, only: test_parameters_t
  use unit_system_oct_m,     only: unit_one, unit_angstrom
 
  use kmeans_clustering_oct_m
 
  implicit none
  private
  public :: test_weighted_kmeans
 
contains
 
  subroutine write_xyz(namespace, fname,  position, species)
    type(namespace_t),       intent(in) :: namespace
    character(len=*),        intent(in) :: fname
    real(real64),            intent(in) :: position(:, :)
    character(len=*),        intent(in) :: species(:)
 
    integer :: iunit, iatom
 
    iunit = io_open(trim(fname)//'.xyz', namespace, action='write', position='asis')
 
    if (mpi_world%is_root()) then
      write(iunit, '(i6)') size(position, 2)
      write(iunit, *)
 
      do iatom = 1, size(position, 2)
        write(iunit, '(A, 1X, 3(F11.6, 1X))') trim(species(iatom)), position(:, iatom) / unit_angstrom%factor
      end do
    endif
 
    call io_close(iunit)
 
  end subroutine write_xyz
 
  subroutine test_weighted_kmeans(namespace)
    type(namespace_t),       intent(in) :: namespace
 
    type(electrons_t), pointer :: sys
    integer                    :: ierr
 
    ! Kmeans
    integer                       :: n_centroids
    real(real64)                  :: inertia
    integer(int64),   allocatable :: seeds(:)
    real(real64),     allocatable :: centroids(:, :)
    character(len=1), allocatable :: dummy_species(:)
 
    push_sub(test_weighted_kmeans)
 
    ! Set up a system of electrons
    call calc_mode_par%set_parallelization(p_strategy_states, default=.false.)
    ! If this is false, the code will crash in v_ks_calc
    sys => electrons_t(namespace, generate_epot=.true.)
    call sys%init_parallelization(mpi_world)
    call states_elec_allocate_wfns(sys%st, sys%gr)
 
    ! Guess at density 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 io_mkdir(static_dir)
 
    ! Output the density: reference xc_vxc_inc.F90
    ! Where the second dim of rho is the spin channel
    call dio_function_output(option__outputformat__cube, static_dir, 'rho', namespace, sys%space, sys%gr, &
      sys%st%rho(:, 1), unit_one, ierr, pos=sys%ions%pos, atoms=sys%ions%atom)
 
    ! Test with a small number of centroids
    n_centroids = 20
 
    allocate(dummy_species(n_centroids), source='A')
    ! Random seed fixed for test reproducibility
    allocate(seeds(1), source=int(101, int64))
    safe_allocate(centroids(1:sys%space%dim, 1:n_centroids))
 
    call sample_initial_centroids(sys%gr, centroids, seed_value=seeds(1))
    call write_xyz(namespace, static_dir // "initial_centroids",  centroids, dummy_species)
 
    call weighted_kmeans(sys%space, sys%gr,  sys%st%rho(1:sys%gr%np, 1), centroids, inertia=inertia)
    write(message(1), '(a, f18.12)') "test_weighted_kmeans: inertia", inertia
    call messages_info(1)
 
    dummy_species = 'B'
    call write_xyz(namespace, static_dir // "final_centroids",  centroids, dummy_species)
 
    ! Free data
    deallocate(seeds)
    deallocate(dummy_species)
    safe_deallocate_a(centroids)
    call states_elec_deallocate_wfns(sys%st)
    safe_deallocate_p(sys)
 
    pop_sub(test_weighted_kmeans)
 
  end subroutine test_weighted_kmeans
 
end module test_kmeans_clustering_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: