Octopus
test_kmeans_clustering.F90
Go to the documentation of this file.
1!! Copyright (C) 2024 A. Buccheri
2!!
3!! This program is free software; you can redistribute it and/or modify
4!! it under the terms of the GNU General Public License as published by
5!! the Free Software Foundation; either version 2, or (at your option)
6!! any later version.
7!!
8!! This program is distributed in the hope that it will be useful,
9!! but WITHOUT ANY WARRANTY; without even the implied warranty of
10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11!! GNU General Public License for more details.
12!!
13!! You should have received a copy of the GNU General Public License
14!! along with this program; if not, write to the Free Software
15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
16!! 02110-1301, USA.
17!!
18
19#include "global.h"
20
22 use, intrinsic :: iso_fortran_env
24 use debug_oct_m
26 use global_oct_m
29 use lcao_oct_m, only: lcao_run
31 use mpi_oct_m, only: mpi_world
37
39
40 implicit none
41 private
42 public :: test_weighted_kmeans
43
44contains
45
50 subroutine test_weighted_kmeans(namespace)
51 type(namespace_t), intent(in) :: namespace
52
53 type(electrons_t), pointer :: sys
54 integer :: ierr
55
56 ! Kmeans
57 integer :: n_centroids
58 real(real64) :: inertia
59 integer(int64), allocatable :: seeds(:)
60 real(real64), allocatable :: centroids(:, :)
61 character(len=1), allocatable :: dummy_species(:)
62
63 push_sub(test_weighted_kmeans)
64
65 ! Set up a system of electrons
66 call calc_mode_par%set_parallelization(p_strategy_states, default=.false.)
67 ! If this is false, the code will crash in v_ks_calc
68 sys => electrons_t(namespace, generate_epot=.true.)
69 call sys%init_parallelization(mpi_world)
70 call states_elec_allocate_wfns(sys%st, sys%gr)
71
72 ! Guess at density using lcao
73 call lcao_run(namespace, sys%space, sys%gr, sys%ions, sys%ext_partners, &
74 sys%st, sys%ks, sys%hm, lmm_r = lmm_r_single_atom)
75
77
78 ! Output the density: reference xc_vxc_inc.F90
79 ! Where the second dim of rho is the spin channel
80 call dio_function_output(option__outputformat__cube, static_dir, 'rho', namespace, sys%space, sys%gr, &
81 sys%st%rho(:, 1), unit_one, ierr, pos=sys%ions%pos, atoms=sys%ions%atom)
82
83 ! Test with a small number of centroids
84 n_centroids = 20
85
86 allocate(dummy_species(n_centroids), source='A')
87 ! Random seed fixed for test reproducibility
88 allocate(seeds(1), source=int(101, int64))
89 safe_allocate(centroids(1:sys%space%dim, 1:n_centroids))
90
91 call sample_initial_centroids(sys%gr, centroids, seed_value=seeds(1))
92 call write_standard_xyz_file(namespace, static_dir // "initial_centroids", centroids, dummy_species)
93
94 call weighted_kmeans(sys%space, sys%gr, sys%st%rho(1:sys%gr%np, 1), centroids, inertia=inertia)
95 write(message(1), '(a, f18.12)') "test_weighted_kmeans: inertia", inertia
96 call messages_info(1)
97
98 dummy_species = 'B'
99 call write_standard_xyz_file(namespace, static_dir // "final_centroids", centroids, dummy_species)
100
101 ! Free data
102 deallocate(seeds)
103 deallocate(dummy_species)
104 safe_deallocate_a(centroids)
105 call states_elec_deallocate_wfns(sys%st)
106 safe_deallocate_p(sys)
107
108 pop_sub(test_weighted_kmeans)
109
110 end subroutine test_weighted_kmeans
111
113
114!! Local Variables:
115!! mode: f90
116!! coding: utf-8
117!! End:
This module handles the calculation mode.
type(calc_mode_par_t), public calc_mode_par
Singleton instance of parallel calculation mode.
integer, parameter, public p_strategy_states
parallelization in states
real(real64), parameter, public lmm_r_single_atom
Default local magnetic moments sphere radius for an isolated system.
Definition: global.F90:220
character(len= *), parameter, public static_dir
Definition: global.F90:265
subroutine, public dio_function_output(how, dir, fname, namespace, space, mesh, ff, unit, ierr, pos, atoms, grp, root)
Top-level IO routine for functions defined on the mesh.
subroutine, public write_standard_xyz_file(namespace, fname, pos, species, header)
Write a standard xyz file with atom labels and positions (in Angstrom).
Definition: io.F90:116
subroutine, public io_close(iunit, grp)
Definition: io.F90:466
subroutine, public io_mkdir(fname, namespace, parents)
Definition: io.F90:360
integer function, public io_open(file, namespace, action, status, form, position, die, recl, grp)
Definition: io.F90:401
subroutine, public weighted_kmeans(space, mesh, weight, centroids, n_iter, centroid_tol, discretize, inertia)
Weighted K-means clustering.
subroutine, public sample_initial_centroids(mesh, centroids, seed_value)
Sample initial centroids from the full mesh.
subroutine, public lcao_run(namespace, space, gr, ions, ext_partners, st, ks, hm, st_start, lmm_r)
Definition: lcao.F90:768
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:162
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:600
type(mpi_grp_t), public mpi_world
Definition: mpi.F90:272
subroutine, public states_elec_deallocate_wfns(st)
Deallocates the KS wavefunctions defined within a states_elec_t structure.
subroutine, public states_elec_allocate_wfns(st, mesh, wfs_type, skip, packed)
Allocates the KS wavefunctions defined within a states_elec_t structure.
subroutine, public test_weighted_kmeans(namespace)
Test weighted kmeans algorithm for a finite system.
This module defines the unit system, used for input and output.
type(unit_t), public unit_angstrom
For XYZ files.
type(unit_t), public unit_one
some special units required for particular quantities
Class describing the electron system.
Definition: electrons.F90:220
int true(void)