(git:b279b6b)
qs_fb_filter_matrix_methods.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
9 
10  USE atomic_kind_types, ONLY: atomic_kind_type,&
12  USE dbcsr_api, ONLY: dbcsr_create,&
13  dbcsr_distribution_type,&
14  dbcsr_finalize,&
15  dbcsr_get_info,&
16  dbcsr_get_stored_coordinates,&
17  dbcsr_put_block,&
18  dbcsr_type,&
19  dbcsr_type_no_symmetry
20  USE fermi_utils, ONLY: fermi,&
22  USE kinds, ONLY: default_string_length,&
23  dp,&
24  int_8
25  USE message_passing, ONLY: mp_para_env_type
26  USE particle_types, ONLY: particle_type
30  fb_atomic_halo_list_obj,&
32  fb_atomic_halo_obj,&
39  USE qs_fb_com_tasks_types, ONLY: &
51  fb_matrix_data_obj,&
54  fb_trial_fns_obj
55  USE string_utilities, ONLY: compress,&
56  uppercase
57 #include "./base/base_uses.f90"
58 
59  IMPLICIT NONE
60 
61  PRIVATE
62 
63  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_fb_filter_matrix_methods'
64 
65  PUBLIC :: fb_fltrmat_build, &
67 
68 CONTAINS
69 
70 ! **************************************************************************************************
71 !> \brief Build the filter matrix, with MPI communications happening at each
72 !> step. Less efficient on communication, but more efficient on
73 !> memory usage (compared to fb_fltrmat_build_2)
74 !> \param H_mat : DBCSR system KS matrix
75 !> \param S_mat : DBCSR system overlap matrix
76 !> \param atomic_halos : list of all local atomic halos, each halo gives
77 !> one atomic matrix and contributes to one blk
78 !> col to the filter matrix
79 !> \param trial_fns : the trial functions to be used to shrink the
80 !> size of the new "filtered" basis
81 !> \param para_env : cp2k parallel environment
82 !> \param particle_set : set of all particles in the system
83 !> \param fermi_level : the fermi level used for defining the filter
84 !> function, which is a Fermi-Dirac distribution
85 !> function
86 !> \param filter_temp : the filter temperature used for defining the
87 !> filter function
88 !> \param name : name given to the filter matrix
89 !> \param filter_mat : DBCSR format filter matrix
90 !> \param tolerance : anything less than tolerance is treated as zero
91 !> \author Lianheng Tong (LT) lianheng.tong@kcl.ac.uk
92 ! **************************************************************************************************
93  SUBROUTINE fb_fltrmat_build(H_mat, &
94  S_mat, &
95  atomic_halos, &
96  trial_fns, &
97  para_env, &
98  particle_set, &
99  fermi_level, &
100  filter_temp, &
101  name, &
102  filter_mat, &
103  tolerance)
104  TYPE(dbcsr_type), POINTER :: h_mat, s_mat
105  TYPE(fb_atomic_halo_list_obj), INTENT(IN) :: atomic_halos
106  TYPE(fb_trial_fns_obj), INTENT(IN) :: trial_fns
107  TYPE(mp_para_env_type), POINTER :: para_env
108  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
109  REAL(kind=dp), INTENT(IN) :: fermi_level, filter_temp
110  CHARACTER(LEN=*), INTENT(IN) :: name
111  TYPE(dbcsr_type), POINTER :: filter_mat
112  REAL(kind=dp), INTENT(IN) :: tolerance
113 
114  CHARACTER(LEN=*), PARAMETER :: routinen = 'fb_fltrmat_build'
115 
116  CHARACTER(LEN=32) :: symmetry_string
117  CHARACTER(LEN=default_string_length) :: name_string
118  INTEGER :: handle, iblkcol, ihalo, ikind, &
119  max_nhalos, nblkcols_total, nhalos
120  INTEGER, DIMENSION(:), POINTER :: col_blk_size, dummy_halo_atoms, ntfns, &
121  row_blk_size
122  LOGICAL :: send_data_only
123  TYPE(atomic_kind_type), POINTER :: atomic_kind
124  TYPE(dbcsr_distribution_type) :: dbcsr_dist
125  TYPE(fb_atomic_halo_obj) :: dummy_atomic_halo
126  TYPE(fb_atomic_halo_obj), DIMENSION(:), POINTER :: halos
127 
128  CALL timeset(routinen, handle)
129 
130  NULLIFY (halos, atomic_kind, ntfns, dummy_halo_atoms, row_blk_size, col_blk_size)
131  CALL fb_atomic_halo_nullify(dummy_atomic_halo)
132 
133  ! filter_mat must be of a dissassociated status (i.e. brand new)
134  cpassert(.NOT. ASSOCIATED(filter_mat))
135 
136  ! get trial function information
137  CALL fb_trial_fns_get(trial_fns=trial_fns, &
138  nfunctions=ntfns)
139 
140  ! calculate the row_blk_size and col_blk_size arrays for
141  ! constructing the filter matrix in DBCSR format
142  ! row_blk_size for the filter matrix is the same as H or S
143  CALL dbcsr_get_info(h_mat, &
144  nblkcols_total=nblkcols_total, &
145  row_blk_size=row_blk_size, &
146  distribution=dbcsr_dist)
147  ALLOCATE (col_blk_size(nblkcols_total))
148  col_blk_size = 0
149  DO iblkcol = 1, nblkcols_total
150  atomic_kind => particle_set(iblkcol)%atomic_kind
151  CALL get_atomic_kind(atomic_kind=atomic_kind, &
152  kind_number=ikind)
153  col_blk_size(iblkcol) = ntfns(ikind)
154  END DO
155  ! DO NOT deallocate cbs if gift=.TRUE. as col_blk_sizes will only point to cbs
156  name_string = name
157  CALL compress(name_string)
158  CALL uppercase(name_string)
159  ! the filter matrix is non-square and is always non-symmetric
160  symmetry_string = dbcsr_type_no_symmetry
161  ! create empty filter matrix
162  ALLOCATE (filter_mat)
163  CALL dbcsr_create(matrix=filter_mat, &
164  name=name_string, &
165  dist=dbcsr_dist, &
166  matrix_type=symmetry_string, &
167  row_blk_size=row_blk_size, &
168  col_blk_size=col_blk_size, &
169  nze=0)
170  DEALLOCATE (col_blk_size)
171 
172  CALL fb_atomic_halo_list_get(atomic_halos=atomic_halos, &
173  nhalos=nhalos, &
174  max_nhalos=max_nhalos, &
175  halos=halos)
176 
177  ! create dummy empty atomic halo
178  CALL fb_atomic_halo_create(dummy_atomic_halo)
179  ALLOCATE (dummy_halo_atoms(0))
180  CALL fb_atomic_halo_set(atomic_halo=dummy_atomic_halo, &
181  owner_atom=0, &
182  owner_id_in_halo=0, &
183  natoms=0, &
184  halo_atoms=dummy_halo_atoms, &
185  nelectrons=0, &
186  sorted=.true.)
187 
188  send_data_only = .false.
189 
190  DO ihalo = 1, max_nhalos
191  IF (ihalo > nhalos) THEN
192  send_data_only = .true.
193  END IF
194  ! construct the filter matrix block by block
195  IF (send_data_only) THEN
196  CALL fb_fltrmat_add_blkcol(h_mat, &
197  s_mat, &
198  dummy_atomic_halo, &
199  trial_fns, &
200  para_env, &
201  particle_set, &
202  fermi_level, &
203  filter_temp, &
204  filter_mat, &
205  tolerance)
206  ELSE
207  CALL fb_fltrmat_add_blkcol(h_mat, &
208  s_mat, &
209  halos(ihalo), &
210  trial_fns, &
211  para_env, &
212  particle_set, &
213  fermi_level, &
214  filter_temp, &
215  filter_mat, &
216  tolerance)
217  END IF ! send_data_only
218  END DO
219 
220  ! finalise the filter matrix
221  CALL dbcsr_finalize(filter_mat)
222 
223  ! cleanup
224  CALL fb_atomic_halo_release(dummy_atomic_halo)
225 
226  CALL timestop(handle)
227 
228  END SUBROUTINE fb_fltrmat_build
229 
230 ! **************************************************************************************************
231 !> \brief Build the filter matrix, with MPI communications grouped together.
232 !> More effcient on communication, less efficient on memory (compared
233 !> to fb_fltrmat_build)
234 !> \param H_mat : DBCSR system KS matrix
235 !> \param S_mat : DBCSR system overlap matrix
236 !> \param atomic_halos : list of all local atomic halos, each halo gives
237 !> one atomic matrix and contributes to one blk
238 !> col to the filter matrix
239 !> \param trial_fns : the trial functions to be used to shrink the
240 !> size of the new "filtered" basis
241 !> \param para_env : cp2k parallel environment
242 !> \param particle_set : set of all particles in the system
243 !> \param fermi_level : the fermi level used for defining the filter
244 !> function, which is a Fermi-Dirac distribution
245 !> function
246 !> \param filter_temp : the filter temperature used for defining the
247 !> filter function
248 !> \param name : name given to the filter matrix
249 !> \param filter_mat : DBCSR format filter matrix
250 !> \param tolerance : anything less than tolerance is treated as zero
251 !> \author Lianheng Tong (LT) lianheng.tong@kcl.ac.uk
252 ! **************************************************************************************************
253  SUBROUTINE fb_fltrmat_build_2(H_mat, &
254  S_mat, &
255  atomic_halos, &
256  trial_fns, &
257  para_env, &
258  particle_set, &
259  fermi_level, &
260  filter_temp, &
261  name, &
262  filter_mat, &
263  tolerance)
264  TYPE(dbcsr_type), POINTER :: h_mat, s_mat
265  TYPE(fb_atomic_halo_list_obj), INTENT(IN) :: atomic_halos
266  TYPE(fb_trial_fns_obj), INTENT(IN) :: trial_fns
267  TYPE(mp_para_env_type), POINTER :: para_env
268  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
269  REAL(kind=dp), INTENT(IN) :: fermi_level, filter_temp
270  CHARACTER(LEN=*), INTENT(IN) :: name
271  TYPE(dbcsr_type), POINTER :: filter_mat
272  REAL(kind=dp), INTENT(IN) :: tolerance
273 
274  CHARACTER(LEN=*), PARAMETER :: routinen = 'fb_fltrmat_build_2'
275 
276  CHARACTER(LEN=default_string_length) :: name_string
277  INTEGER :: handle, iblkcol, ihalo, ikind, &
278  natoms_global, natoms_in_halo, &
279  nblkcols_total, nblks_recv, nhalos, &
280  nmax
281  INTEGER, DIMENSION(:), POINTER :: col_blk_size, ntfns, row_blk_size
282  LOGICAL :: check_ok
283  TYPE(atomic_kind_type), POINTER :: atomic_kind
284  TYPE(dbcsr_distribution_type) :: dbcsr_dist
285  TYPE(fb_atomic_halo_obj), DIMENSION(:), POINTER :: halos
286  TYPE(fb_com_atom_pairs_obj) :: atmatrix_blks_recv, atmatrix_blks_send, &
287  filter_mat_blks_recv, &
288  filter_mat_blks_send
289  TYPE(fb_matrix_data_obj) :: filter_mat_data, h_mat_data, s_mat_data
290 
291  CALL timeset(routinen, handle)
292 
293  NULLIFY (halos, atomic_kind, row_blk_size, col_blk_size, ntfns)
294 
295  ! filter_mat must be of a dissassociated status (i.e. brand new)
296  check_ok = .NOT. ASSOCIATED(filter_mat)
297  cpassert(check_ok)
298 
299  ! get total number of atoms
300  natoms_global = SIZE(particle_set)
301 
302  ! get trial function information
303  CALL fb_trial_fns_get(trial_fns=trial_fns, &
304  nfunctions=ntfns)
305 
306  ! calculate the row_blk_size and col_blk_size arrays for
307  ! constructing the filter matrix in DBCSR format
308  ! row_blk_size for the filter matrix is the same as H or S
309  CALL dbcsr_get_info(h_mat, &
310  nblkcols_total=nblkcols_total, &
311  row_blk_size=row_blk_size, &
312  distribution=dbcsr_dist)
313  ALLOCATE (col_blk_size(nblkcols_total))
314  col_blk_size = 0
315  DO iblkcol = 1, nblkcols_total
316  atomic_kind => particle_set(iblkcol)%atomic_kind
317  CALL get_atomic_kind(atomic_kind=atomic_kind, &
318  kind_number=ikind)
319  col_blk_size(iblkcol) = ntfns(ikind)
320  END DO
321  ! DO NOT deallocate cbs if gift=.TRUE. as col_blk_sizes will only point to cbs
322  name_string = name
323  CALL compress(name_string)
324  CALL uppercase(name_string)
325  ! create empty filter matrix (it is always non-symmetric as it is non-square)
326  ALLOCATE (filter_mat)
327  CALL dbcsr_create(matrix=filter_mat, &
328  name=name_string, &
329  dist=dbcsr_dist, &
330  matrix_type=dbcsr_type_no_symmetry, &
331  row_blk_size=row_blk_size, &
332  col_blk_size=col_blk_size, &
333  nze=0)
334  DEALLOCATE (col_blk_size)
335 
336  ! get all the blocks required for constructing atomic matrics, and
337  ! store it in a fb_matrix_data object
338  CALL fb_matrix_data_nullify(h_mat_data)
339  CALL fb_matrix_data_nullify(s_mat_data)
340  CALL fb_com_atom_pairs_nullify(atmatrix_blks_send)
341  CALL fb_com_atom_pairs_nullify(atmatrix_blks_recv)
342  CALL fb_com_atom_pairs_create(atmatrix_blks_send)
343  CALL fb_com_atom_pairs_create(atmatrix_blks_recv)
344  ! H matrix
346  atomic_halos, &
347  para_env, &
348  atmatrix_blks_send, &
349  atmatrix_blks_recv)
350  CALL fb_com_atom_pairs_get(atom_pairs=atmatrix_blks_recv, &
351  npairs=nblks_recv)
352  CALL fb_matrix_data_create(h_mat_data, &
353  nblks_recv, &
354  natoms_global)
355  CALL fb_com_atom_pairs_gather_blks(h_mat, &
356  atmatrix_blks_send, &
357  atmatrix_blks_recv, &
358  para_env, &
359  h_mat_data)
360  ! S matrix
362  atomic_halos, &
363  para_env, &
364  atmatrix_blks_send, &
365  atmatrix_blks_recv)
366  CALL fb_com_atom_pairs_get(atom_pairs=atmatrix_blks_recv, &
367  npairs=nblks_recv)
368  CALL fb_matrix_data_create(s_mat_data, &
369  nblks_recv, &
370  natoms_global)
371  CALL fb_com_atom_pairs_gather_blks(s_mat, &
372  atmatrix_blks_send, &
373  atmatrix_blks_recv, &
374  para_env, &
375  s_mat_data)
376  ! cleanup
377  CALL fb_com_atom_pairs_release(atmatrix_blks_send)
378  CALL fb_com_atom_pairs_release(atmatrix_blks_recv)
379 
380  ! make filter matrix blocks one by one and store in an
381  ! matrix_data_obj
382  CALL fb_matrix_data_nullify(filter_mat_data)
383  CALL fb_atomic_halo_list_get(atomic_halos=atomic_halos, &
384  nhalos=nhalos, &
385  halos=halos)
386  nmax = 0
387  DO ihalo = 1, nhalos
388  CALL fb_atomic_halo_get(atomic_halo=halos(ihalo), &
389  natoms=natoms_in_halo)
390  nmax = nmax + natoms_in_halo
391  END DO
392  CALL fb_matrix_data_create(filter_mat_data, &
393  nmax, &
394  natoms_global)
395  DO ihalo = 1, nhalos
396  CALL fb_fltrmat_add_blkcol_2(h_mat, &
397  s_mat, &
398  h_mat_data, &
399  s_mat_data, &
400  halos(ihalo), &
401  trial_fns, &
402  particle_set, &
403  fermi_level, &
404  filter_temp, &
405  filter_mat_data, &
406  tolerance)
407  END DO
408  ! clean up
409  CALL fb_matrix_data_release(h_mat_data)
410  CALL fb_matrix_data_release(s_mat_data)
411 
412  ! distribute the relevant blocks from the matrix_data_obj to DBCSR
413  ! filter matrix
414  CALL fb_com_atom_pairs_nullify(filter_mat_blks_send)
415  CALL fb_com_atom_pairs_nullify(filter_mat_blks_recv)
416  CALL fb_com_atom_pairs_create(filter_mat_blks_send)
417  CALL fb_com_atom_pairs_create(filter_mat_blks_recv)
418  CALL fb_fltrmat_generate_com_pairs_2(filter_mat, &
419  atomic_halos, &
420  para_env, &
421  filter_mat_blks_send, &
422  filter_mat_blks_recv)
423  CALL fb_com_atom_pairs_distribute_blks(filter_mat_data, &
424  filter_mat_blks_send, &
425  filter_mat_blks_recv, &
426  para_env, &
427  filter_mat)
428  ! cleanup
429  CALL fb_com_atom_pairs_release(filter_mat_blks_send)
430  CALL fb_com_atom_pairs_release(filter_mat_blks_recv)
431  CALL fb_matrix_data_release(filter_mat_data)
432 
433  ! finalise matrix
434  CALL dbcsr_finalize(filter_mat)
435 
436  CALL timestop(handle)
437 
438  END SUBROUTINE fb_fltrmat_build_2
439 
440 ! **************************************************************************************************
441 !> \brief Add a computed blocks in one column to the filter matrix. This
442 !> version is used by fb_fltrmat_build, for the case where MPI
443 !> communications are done at each step
444 !> It does not finalise the filter matrix
445 !> \param H_mat : DBCSR system KS matrix
446 !> \param S_mat : DBCSR system overlap matrix
447 !> \param atomic_halo : the halo that contributes to the blk
448 !> col of the filter matrix
449 !> \param trial_fns ...
450 !> \param para_env : cp2k parallel environment
451 !> \param particle_set : set of all particles in the system
452 !> \param fermi_level : the fermi level used for defining the filter
453 !> function, which is a Fermi-Dirac distribution
454 !> function
455 !> \param filter_temp : the filter temperature used for defining the
456 !> filter function
457 !> \param filter_mat : DBCSR format filter matrix
458 !> \param tolerance : anything smaller than tolerance is treated as zero
459 !> \author Lianheng Tong (LT) lianheng.tong@kcl.ac.uk
460 ! **************************************************************************************************
461  SUBROUTINE fb_fltrmat_add_blkcol(H_mat, &
462  S_mat, &
463  atomic_halo, &
464  trial_fns, &
465  para_env, &
466  particle_set, &
467  fermi_level, &
468  filter_temp, &
469  filter_mat, &
470  tolerance)
471  TYPE(dbcsr_type), POINTER :: h_mat, s_mat
472  TYPE(fb_atomic_halo_obj), INTENT(IN) :: atomic_halo
473  TYPE(fb_trial_fns_obj), INTENT(IN) :: trial_fns
474  TYPE(mp_para_env_type), POINTER :: para_env
475  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
476  REAL(kind=dp), INTENT(IN) :: fermi_level, filter_temp
477  TYPE(dbcsr_type), POINTER :: filter_mat
478  REAL(kind=dp), INTENT(IN) :: tolerance
479 
480  CHARACTER(LEN=*), PARAMETER :: routinen = 'fb_fltrmat_add_blkcol'
481 
482  INTEGER :: handle, handle_mpi, iatom_global, iatom_in_halo, ind, ipair, ipe, itrial, &
483  jatom_global, jatom_in_halo, jkind, natoms_global, natoms_in_halo, ncols_atmatrix, &
484  ncols_blk, nrows_atmatrix, nrows_blk, numprocs, pe, recv_encode, send_encode, stat
485  INTEGER(KIND=int_8), DIMENSION(:), POINTER :: pairs_recv, pairs_send
486  INTEGER, ALLOCATABLE, DIMENSION(:) :: atomic_h_blk_col_start, atomic_h_blk_row_start, &
487  atomic_s_blk_col_start, atomic_s_blk_row_start, col_block_size_data, ind_in_halo, &
488  recv_disps, recv_pair_count, recv_pair_disps, recv_sizes, send_disps, send_pair_count, &
489  send_pair_disps, send_sizes
490  INTEGER, DIMENSION(:), POINTER :: halo_atoms, ntfns, row_block_size_data
491  INTEGER, DIMENSION(:, :), POINTER :: tfns
492  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: recv_buf, send_buf
493  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: atomic_filter_mat, atomic_h, atomic_s
494  TYPE(atomic_kind_type), POINTER :: atomic_kind
495  TYPE(fb_com_atom_pairs_obj) :: com_pairs_recv, com_pairs_send
496 
497  CALL timeset(routinen, handle)
498 
499  NULLIFY (atomic_kind, halo_atoms, ntfns, pairs_send, pairs_recv, &
500  row_block_size_data, tfns)
501  CALL fb_com_atom_pairs_nullify(com_pairs_send)
502  CALL fb_com_atom_pairs_nullify(com_pairs_recv)
503 
504  ! ----------------------------------------------------------------------
505  ! Get communication buffers ready
506  ! ----------------------------------------------------------------------
507 
508  ! generate send and recv atom pairs
509  CALL fb_com_atom_pairs_create(com_pairs_send)
510  CALL fb_com_atom_pairs_create(com_pairs_recv)
511  CALL fb_fltrmat_generate_com_pairs(filter_mat, &
512  atomic_halo, &
513  para_env, &
514  com_pairs_send, &
515  com_pairs_recv)
516  CALL fb_com_atom_pairs_get(atom_pairs=com_pairs_send, &
517  natoms_encode=send_encode, &
518  pairs=pairs_send)
519  CALL fb_com_atom_pairs_get(atom_pairs=com_pairs_recv, &
520  natoms_encode=recv_encode, &
521  pairs=pairs_recv)
522 
523  ! get para_env info
524  numprocs = para_env%num_pe
525  ! me = para_env%mepos + 1 ! my process id, starting counting from 1
526 
527  ! obtain trail function information
528  CALL fb_trial_fns_get(trial_fns=trial_fns, &
529  nfunctions=ntfns, &
530  functions=tfns)
531 
532  ! obtain row and col block size data for filter matrix
533  CALL dbcsr_get_info(h_mat, row_blk_size=row_block_size_data)
534  natoms_global = SIZE(particle_set)
535  ALLOCATE (col_block_size_data(natoms_global))
536  DO jatom_global = 1, natoms_global
537  atomic_kind => particle_set(jatom_global)%atomic_kind
538  CALL get_atomic_kind(atomic_kind=atomic_kind, kind_number=jkind)
539  col_block_size_data(jatom_global) = ntfns(jkind)
540  END DO
541 
542  ! allocate temporary arrays for send
543  ALLOCATE (send_sizes(numprocs))
544  ALLOCATE (send_disps(numprocs))
545  ALLOCATE (send_pair_count(numprocs))
546  ALLOCATE (send_pair_disps(numprocs))
547  ! setup send buffer sizes
548  CALL fb_com_atom_pairs_calc_buffer_sizes(com_pairs_send, &
549  numprocs, &
550  row_block_size_data, &
551  col_block_size_data, &
552  send_sizes, &
553  send_disps, &
554  send_pair_count, &
555  send_pair_disps)
556  ! allocate send buffer
557  ALLOCATE (send_buf(sum(send_sizes)))
558 
559  ! allocate temporary array for recv
560  ALLOCATE (recv_sizes(numprocs))
561  ALLOCATE (recv_disps(numprocs))
562  ALLOCATE (recv_pair_count(numprocs))
563  ALLOCATE (recv_pair_disps(numprocs))
564  ! setup recv buffer sizes
565  CALL fb_com_atom_pairs_calc_buffer_sizes(com_pairs_recv, &
566  numprocs, &
567  row_block_size_data, &
568  col_block_size_data, &
569  recv_sizes, &
570  recv_disps, &
571  recv_pair_count, &
572  recv_pair_disps)
573  ! allocate recv buffer
574  ALLOCATE (recv_buf(sum(recv_sizes)))
575 
576  ! ----------------------------------------------------------------------
577  ! Construct atomic filter matrix for this atomic_halo
578  ! ----------------------------------------------------------------------
579 
580  CALL fb_atomic_halo_get(atomic_halo=atomic_halo, &
581  natoms=natoms_in_halo, &
582  halo_atoms=halo_atoms)
583 
584  ! construct atomic matrix for H for atomic_halo
585  ALLOCATE (atomic_h_blk_row_start(natoms_in_halo + 1), &
586  atomic_h_blk_col_start(natoms_in_halo + 1), &
587  stat=stat)
588  cpassert(stat == 0)
589  CALL fb_atmatrix_calc_size(h_mat, &
590  atomic_halo, &
591  nrows_atmatrix, &
592  ncols_atmatrix, &
593  atomic_h_blk_row_start, &
594  atomic_h_blk_col_start)
595 
596  ALLOCATE (atomic_h(nrows_atmatrix, ncols_atmatrix))
597  CALL fb_atmatrix_construct(h_mat, &
598  atomic_halo, &
599  para_env, &
600  atomic_h, &
601  atomic_h_blk_row_start, &
602  atomic_h_blk_col_start)
603 
604  ! construct atomic matrix for S for atomic_halo
605  ALLOCATE (atomic_s_blk_row_start(natoms_in_halo + 1), &
606  atomic_s_blk_col_start(natoms_in_halo + 1), &
607  stat=stat)
608  cpassert(stat == 0)
609  CALL fb_atmatrix_calc_size(s_mat, &
610  atomic_halo, &
611  nrows_atmatrix, &
612  ncols_atmatrix, &
613  atomic_s_blk_row_start, &
614  atomic_s_blk_col_start)
615  ALLOCATE (atomic_s(nrows_atmatrix, ncols_atmatrix))
616  CALL fb_atmatrix_construct(s_mat, &
617  atomic_halo, &
618  para_env, &
619  atomic_s, &
620  atomic_s_blk_row_start, &
621  atomic_s_blk_col_start)
622 
623  ! construct the atomic filter matrix
624  ALLOCATE (atomic_filter_mat(nrows_atmatrix, ncols_atmatrix))
625  ! calculate atomic filter matrix only if it is non-zero sized
626  IF (nrows_atmatrix > 0 .AND. ncols_atmatrix > 0) THEN
627  CALL fb_fltrmat_build_atomic_fltrmat(atomic_h, &
628  atomic_s, &
629  fermi_level, &
630  filter_temp, &
631  atomic_filter_mat, &
632  tolerance)
633  END IF
634 
635  ! ----------------------------------------------------------------------
636  ! Construct filter matrix blocks and add to the correct locations
637  ! in send_buffer
638  ! ----------------------------------------------------------------------
639 
640  ! preconstruct iatom_global to iatom_in_halo map
641  ALLOCATE (ind_in_halo(natoms_global))
642  ind_in_halo = 0
643  DO iatom_in_halo = 1, natoms_in_halo
644  iatom_global = halo_atoms(iatom_in_halo)
645  ind_in_halo(iatom_global) = iatom_in_halo
646  END DO
647 
648  ! initialise send buffer
649  IF (SIZE(send_buf) > 0) send_buf = 0.0_dp
650  ! assign values
651  DO ipe = 1, numprocs
652  send_sizes(ipe) = 0
653  DO ipair = 1, send_pair_count(ipe)
654  CALL fb_com_atom_pairs_decode(pairs_send(send_pair_disps(ipe) + ipair), &
655  pe, iatom_global, jatom_global, &
656  send_encode)
657  iatom_in_halo = ind_in_halo(iatom_global)
658  cpassert(iatom_in_halo > 0)
659  jatom_in_halo = ind_in_halo(jatom_global)
660  cpassert(jatom_in_halo > 0)
661  atomic_kind => particle_set(jatom_global)%atomic_kind
662  CALL get_atomic_kind(atomic_kind=atomic_kind, &
663  kind_number=jkind)
664  nrows_blk = row_block_size_data(iatom_global)
665  ncols_blk = ntfns(jkind)
666 
667  ! do it column-wise one trial function at a time
668  DO itrial = 1, ntfns(jkind)
669  ind = send_disps(ipe) + send_sizes(ipe) + (itrial - 1)*nrows_blk
670  CALL dgemv("N", &
671  nrows_blk, &
672  ncols_atmatrix, &
673  1.0_dp, &
674  atomic_filter_mat( &
675  atomic_h_blk_row_start(iatom_in_halo): &
676  atomic_h_blk_row_start(iatom_in_halo + 1) - 1, &
677  1:ncols_atmatrix &
678  ), &
679  nrows_blk, &
680  atomic_s( &
681  1:nrows_atmatrix, &
682  atomic_s_blk_col_start(jatom_in_halo) + &
683  tfns(itrial, jkind) - 1 &
684  ), &
685  1, &
686  0.0_dp, &
687  send_buf(ind + 1:ind + nrows_blk), &
688  1)
689  END DO ! itrial
690  send_sizes(ipe) = send_sizes(ipe) + nrows_blk*ncols_blk
691  END DO ! ipair
692  END DO ! ipe
693 
694  DEALLOCATE (atomic_h)
695  DEALLOCATE (atomic_h_blk_row_start)
696  DEALLOCATE (atomic_s)
697  DEALLOCATE (atomic_s_blk_row_start)
698  DEALLOCATE (atomic_filter_mat)
699  DEALLOCATE (ind_in_halo)
700 
701  ! ----------------------------------------------------------------------
702  ! Do communication
703  ! ----------------------------------------------------------------------
704 
705  CALL timeset("fb_fltrmat_add_blkcol_mpi", handle_mpi)
706 
707  CALL para_env%alltoall(send_buf, send_sizes, send_disps, &
708  recv_buf, recv_sizes, recv_disps)
709 
710  CALL timestop(handle_mpi)
711 
712  DEALLOCATE (send_buf)
713  DEALLOCATE (send_sizes)
714  DEALLOCATE (send_disps)
715  DEALLOCATE (send_pair_count)
716  DEALLOCATE (send_pair_disps)
717 
718  ! ----------------------------------------------------------------------
719  ! Unpack the recv buffer and add the blocks to correct parts of
720  ! the DBCSR filter matrix
721  ! ----------------------------------------------------------------------
722 
723  DO ipe = 1, numprocs
724  recv_sizes(ipe) = 0
725  DO ipair = 1, recv_pair_count(ipe)
726  CALL fb_com_atom_pairs_decode(pairs_recv(recv_pair_disps(ipe) + ipair), &
727  pe, iatom_global, jatom_global, &
728  recv_encode)
729  nrows_blk = row_block_size_data(iatom_global)
730  ncols_blk = col_block_size_data(jatom_global)
731  ind = recv_disps(ipe) + recv_sizes(ipe)
732  CALL dbcsr_put_block(filter_mat, &
733  iatom_global, jatom_global, &
734  recv_buf((ind + 1):(ind + nrows_blk*ncols_blk)))
735  recv_sizes(ipe) = recv_sizes(ipe) + nrows_blk*ncols_blk
736  END DO ! ipair
737  END DO ! ipe
738 
739  ! cleanup rest of the temporary arrays
740  DEALLOCATE (recv_buf)
741  DEALLOCATE (recv_sizes)
742  DEALLOCATE (recv_pair_count)
743  DEALLOCATE (recv_pair_disps)
744 
745  CALL fb_com_atom_pairs_release(com_pairs_send)
746  CALL fb_com_atom_pairs_release(com_pairs_recv)
747 
748  ! cannot finalise the matrix until all blocks has been added
749 
750  CALL timestop(handle)
751 
752  END SUBROUTINE fb_fltrmat_add_blkcol
753 
754 ! **************************************************************************************************
755 !> \brief Computed blocks in one filter matrix column. This version is used by
756 !> fb_fltrmat_build_2, where MPI communication is done collectively
757 !> \param H_mat : DBCSR system KS matrix
758 !> \param S_mat : DBCSR system overlap matrix
759 !> \param H_mat_data : local storage of the relevant H_mat matrix blocks
760 !> \param S_mat_data : local storage of the relevant S_mat matrix blocks
761 !> \param atomic_halo : the halo that contributes to the blk
762 !> col of the filter matrix
763 !> \param trial_fns : trial functions data
764 !> \param particle_set : set of all particles in the system
765 !> \param fermi_level : the fermi level used for defining the filter
766 !> function, which is a Fermi-Dirac distribution
767 !> function
768 !> \param filter_temp : the filter temperature used for defining the
769 !> filter function
770 !> \param filter_mat_data : local storage for the the computed filter matrix
771 !> blocks
772 !> \param tolerance : anything less than this is regarded as zero
773 !> \author Lianheng Tong (LT) lianheng.tong@kcl.ac.uk
774 ! **************************************************************************************************
775  SUBROUTINE fb_fltrmat_add_blkcol_2(H_mat, &
776  S_mat, &
777  H_mat_data, &
778  S_mat_data, &
779  atomic_halo, &
780  trial_fns, &
781  particle_set, &
782  fermi_level, &
783  filter_temp, &
784  filter_mat_data, &
785  tolerance)
786  TYPE(dbcsr_type), POINTER :: h_mat, s_mat
787  TYPE(fb_matrix_data_obj), INTENT(IN) :: h_mat_data, s_mat_data
788  TYPE(fb_atomic_halo_obj), INTENT(IN) :: atomic_halo
789  TYPE(fb_trial_fns_obj), INTENT(IN) :: trial_fns
790  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
791  REAL(kind=dp), INTENT(IN) :: fermi_level, filter_temp
792  TYPE(fb_matrix_data_obj), INTENT(INOUT) :: filter_mat_data
793  REAL(kind=dp), INTENT(IN) :: tolerance
794 
795  CHARACTER(LEN=*), PARAMETER :: routinen = 'fb_fltrmat_add_blkcol_2'
796 
797  INTEGER :: handle, iatom_global, iatom_in_halo, itrial, jatom_global, jatom_in_halo, jkind, &
798  natoms_global, natoms_in_halo, ncols_atmatrix, ncols_blk, ncols_blk_max, nrows_atmatrix, &
799  nrows_blk, nrows_blk_max, stat
800  INTEGER, ALLOCATABLE, DIMENSION(:) :: atomic_h_blk_col_start, atomic_h_blk_row_start, &
801  atomic_s_blk_col_start, atomic_s_blk_row_start, col_block_size_data
802  INTEGER, DIMENSION(:), POINTER :: halo_atoms, ntfns, row_block_size_data
803  INTEGER, DIMENSION(:, :), POINTER :: tfns
804  LOGICAL :: check_ok
805  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: atomic_filter_mat, atomic_h, atomic_s, &
806  mat_blk
807  TYPE(atomic_kind_type), POINTER :: atomic_kind
808 
809  CALL timeset(routinen, handle)
810 
811  NULLIFY (atomic_kind, halo_atoms, ntfns, row_block_size_data, tfns)
812 
813  check_ok = fb_matrix_data_has_data(h_mat_data)
814  cpassert(check_ok)
815  check_ok = fb_matrix_data_has_data(s_mat_data)
816  cpassert(check_ok)
817 
818  ! obtain trial function information
819  CALL fb_trial_fns_get(trial_fns=trial_fns, &
820  nfunctions=ntfns, &
821  functions=tfns)
822 
823  ! obtain row and col block size data for filter matrix
824  CALL dbcsr_get_info(h_mat, row_blk_size=row_block_size_data)
825  natoms_global = SIZE(particle_set)
826  ALLOCATE (col_block_size_data(natoms_global))
827  DO jatom_global = 1, natoms_global
828  atomic_kind => particle_set(jatom_global)%atomic_kind
829  CALL get_atomic_kind(atomic_kind=atomic_kind, kind_number=jkind)
830  col_block_size_data(jatom_global) = ntfns(jkind)
831  END DO
832 
833  ! ----------------------------------------------------------------------
834  ! Construct atomic filter matrix for this atomic_halo
835  ! ----------------------------------------------------------------------
836 
837  CALL fb_atomic_halo_get(atomic_halo=atomic_halo, &
838  natoms=natoms_in_halo, &
839  halo_atoms=halo_atoms)
840 
841  ! construct atomic matrix for H for atomic_halo
842  ALLOCATE (atomic_h_blk_row_start(natoms_in_halo + 1), &
843  atomic_h_blk_col_start(natoms_in_halo + 1), &
844  stat=stat)
845  cpassert(stat == 0)
846  CALL fb_atmatrix_calc_size(h_mat, &
847  atomic_halo, &
848  nrows_atmatrix, &
849  ncols_atmatrix, &
850  atomic_h_blk_row_start, &
851  atomic_h_blk_col_start)
852  ALLOCATE (atomic_h(nrows_atmatrix, ncols_atmatrix))
853  CALL fb_atmatrix_construct_2(h_mat_data, &
854  atomic_halo, &
855  atomic_h, &
856  atomic_h_blk_row_start, &
857  atomic_h_blk_col_start)
858 
859  ! construct atomic matrix for S for atomic_halo
860  ALLOCATE (atomic_s_blk_row_start(natoms_in_halo + 1), &
861  atomic_s_blk_col_start(natoms_in_halo + 1), &
862  stat=stat)
863  cpassert(stat == 0)
864  CALL fb_atmatrix_calc_size(s_mat, &
865  atomic_halo, &
866  nrows_atmatrix, &
867  ncols_atmatrix, &
868  atomic_s_blk_row_start, &
869  atomic_s_blk_col_start)
870  ALLOCATE (atomic_s(nrows_atmatrix, ncols_atmatrix))
871  CALL fb_atmatrix_construct_2(s_mat_data, &
872  atomic_halo, &
873  atomic_s, &
874  atomic_s_blk_row_start, &
875  atomic_s_blk_col_start)
876 
877  ! construct the atomic filter matrix
878  ALLOCATE (atomic_filter_mat(nrows_atmatrix, ncols_atmatrix))
879  ! calculate atomic filter matrix only if it is non-zero sized
880  IF (nrows_atmatrix > 0 .AND. ncols_atmatrix > 0) THEN
881  CALL fb_fltrmat_build_atomic_fltrmat(atomic_h, &
882  atomic_s, &
883  fermi_level, &
884  filter_temp, &
885  atomic_filter_mat, &
886  tolerance)
887  END IF
888 
889  ! ----------------------------------------------------------------------
890  ! Construct filter matrix block and add to filter_mat_data
891  ! ----------------------------------------------------------------------
892 
893  CALL fb_atomic_halo_get(atomic_halo=atomic_halo, &
894  owner_atom=jatom_global, &
895  owner_id_in_halo=jatom_in_halo)
896  nrows_blk_max = maxval(row_block_size_data)
897  ncols_blk_max = maxval(ntfns)
898  ALLOCATE (mat_blk(nrows_blk_max, ncols_blk_max))
899  mat_blk(:, :) = 0.0_dp
900  DO iatom_in_halo = 1, natoms_in_halo
901  iatom_global = halo_atoms(iatom_in_halo)
902  atomic_kind => particle_set(jatom_global)%atomic_kind
903  CALL get_atomic_kind(atomic_kind=atomic_kind, &
904  kind_number=jkind)
905  nrows_blk = row_block_size_data(iatom_global)
906  ncols_blk = ntfns(jkind)
907 
908  ! ALLOCATE(mat_blk(nrows_blk,ncols_blk) STAT=stat)
909  ! CPPostcondition(stat==0, cp_failure_level, routineP,failure)
910 
911  ! do it column-wise one trial function at a time
912  DO itrial = 1, ntfns(jkind)
913  CALL dgemv("N", &
914  nrows_blk, &
915  ncols_atmatrix, &
916  1.0_dp, &
917  atomic_filter_mat( &
918  atomic_h_blk_row_start(iatom_in_halo): &
919  atomic_h_blk_row_start(iatom_in_halo + 1) - 1, &
920  1:ncols_atmatrix &
921  ), &
922  nrows_blk, &
923  atomic_s( &
924  1:nrows_atmatrix, &
925  atomic_s_blk_col_start(jatom_in_halo) + &
926  tfns(itrial, jkind) - 1 &
927  ), &
928  1, &
929  0.0_dp, &
930  mat_blk( &
931  1:nrows_blk, &
932  itrial), &
933  1)
934  END DO ! itrial
935  CALL fb_matrix_data_add(filter_mat_data, &
936  iatom_global, &
937  jatom_global, &
938  mat_blk(1:nrows_blk, 1:ncols_blk))
939  ! DEALLOCATE(mat_blk, STAT=stat)
940  ! CPPostcondition(stat==0, cp_failure_level, routineP,failure)
941  END DO ! iatom_in_halo
942  DEALLOCATE (mat_blk)
943 
944  ! clean up
945  DEALLOCATE (atomic_h)
946  DEALLOCATE (atomic_h_blk_row_start)
947  DEALLOCATE (atomic_s)
948  DEALLOCATE (atomic_s_blk_row_start)
949  DEALLOCATE (atomic_filter_mat)
950 
951  CALL timestop(handle)
952 
953  END SUBROUTINE fb_fltrmat_add_blkcol_2
954 
955 ! **************************************************************************************************
956 !> \brief generate the list of blocks (atom pairs) to be sent and received
957 !> in order to construct the filter matrix for each atomic halo.
958 !> This version is for use with fb_fltrmat_build, where MPI
959 !> communications are done at each step
960 !> \param filter_mat : DBCSR formatted filter matrix
961 !> \param atomic_halo : the halo that contributes to a blk
962 !> col of the filter matrix
963 !> \param para_env : cp2k parallel environment
964 !> \param atom_pairs_send : list of blocks to be sent
965 !> \param atom_pairs_recv : list of blocks to be received
966 !> \author Lianheng Tong (LT) lianheng.tong@kcl.ac.uk
967 ! **************************************************************************************************
968  SUBROUTINE fb_fltrmat_generate_com_pairs(filter_mat, &
969  atomic_halo, &
970  para_env, &
971  atom_pairs_send, &
972  atom_pairs_recv)
973  TYPE(dbcsr_type), POINTER :: filter_mat
974  TYPE(fb_atomic_halo_obj), INTENT(IN) :: atomic_halo
975  TYPE(mp_para_env_type), POINTER :: para_env
976  TYPE(fb_com_atom_pairs_obj), INTENT(INOUT) :: atom_pairs_send, atom_pairs_recv
977 
978  CHARACTER(LEN=*), PARAMETER :: routinen = 'fb_fltrmat_generate_com_pairs'
979 
980  INTEGER :: dest, handle, iatom_global, &
981  iatom_in_halo, itask, jatom_global, &
982  natoms_in_halo, nblkrows_total, &
983  ntasks_send
984  INTEGER(KIND=int_8), DIMENSION(:, :), POINTER :: tasks_send
985  INTEGER, DIMENSION(:), POINTER :: halo_atoms
986  TYPE(fb_com_tasks_obj) :: com_tasks_recv, com_tasks_send
987 
988  CALL timeset(routinen, handle)
989 
990  NULLIFY (tasks_send)
991  CALL fb_com_tasks_nullify(com_tasks_send)
992  CALL fb_com_tasks_nullify(com_tasks_recv)
993 
994  ! initialise atom_pairs_send and atom_pairs_recv
995  IF (fb_com_atom_pairs_has_data(atom_pairs_send)) THEN
996  CALL fb_com_atom_pairs_init(atom_pairs_send)
997  ELSE
998  CALL fb_com_atom_pairs_create(atom_pairs_send)
999  END IF
1000  IF (fb_com_atom_pairs_has_data(atom_pairs_recv)) THEN
1001  CALL fb_com_atom_pairs_init(atom_pairs_recv)
1002  ELSE
1003  CALL fb_com_atom_pairs_create(atom_pairs_recv)
1004  END IF
1005 
1006  ! The total number of filter matrix blocks each processor is going
1007  ! to construct equals to the total number of halo atoms in all of
1008  ! the atomic halos local to the processor. The number of send
1009  ! tasks will not exceed this. We do one halo (col) at a time, and
1010  ! each call of this subroutine will only work on one filter matrix
1011  ! col corresponding to atomic_halo.
1012 
1013  ! The col atom block index for each filter matrix block are the
1014  ! owner atom of each halo. The row atom block index for each
1015  ! filter matrix block corresponding to each col are the halo atoms
1016  ! of the corresponding halos. Filter matrix is non-symmetric: it
1017  ! is non-square, because the blocks themselves are non-square
1018 
1019  CALL fb_atomic_halo_get(atomic_halo=atomic_halo, &
1020  owner_atom=jatom_global, &
1021  natoms=natoms_in_halo, &
1022  halo_atoms=halo_atoms)
1023  ntasks_send = natoms_in_halo
1024 
1025  ! allocate send tasks
1026  ALLOCATE (tasks_send(task_n_records, ntasks_send))
1027 
1028  ! Get the total number of atoms, this can be obtained from the
1029  ! total number of block rows in the DBCSR filter matrix. We
1030  ! assumes that before calling this subroutine, the filter_mat has
1031  ! already been created and initialised: i.e. using
1032  ! dbcsr_create_new. Even if the matrix is at the moment empty,
1033  ! the attribute nblkrows_total is already assigned from the dbcsr
1034  ! distribution data
1035  CALL dbcsr_get_info(filter_mat, &
1036  nblkrows_total=nblkrows_total)
1037 
1038  ! source is always the local processor
1039  associate(src => para_env%mepos)
1040  ! construct send tasks
1041  itask = 1
1042  DO iatom_in_halo = 1, natoms_in_halo
1043  iatom_global = halo_atoms(iatom_in_halo)
1044  ! find where the constructed block of filter matrix belongs to
1045  CALL dbcsr_get_stored_coordinates(filter_mat, &
1046  iatom_global, &
1047  jatom_global, &
1048  processor=dest)
1049  ! create the send tasks
1050  tasks_send(task_dest, itask) = dest
1051  tasks_send(task_src, itask) = src
1052  CALL fb_com_tasks_encode_pair(tasks_send(task_pair, itask), &
1053  iatom_global, jatom_global, &
1054  nblkrows_total)
1055  ! calculation of cost not implemented at the moment
1056  tasks_send(task_cost, itask) = 0
1057  itask = itask + 1
1058  END DO ! iatom_in_halo
1059  END associate
1060 
1061  CALL fb_com_tasks_create(com_tasks_recv)
1062  CALL fb_com_tasks_create(com_tasks_send)
1063 
1064  CALL fb_com_tasks_set(com_tasks=com_tasks_send, &
1065  task_dim=task_n_records, &
1066  ntasks=ntasks_send, &
1067  nencode=nblkrows_total, &
1068  tasks=tasks_send)
1069 
1070  ! generate the recv task list (tasks_recv) from the send task list
1071  CALL fb_com_tasks_transpose_dest_src(com_tasks_recv, "<", com_tasks_send, &
1072  para_env)
1073 
1074  ! task lists are now complete, now construct the atom_pairs_send
1075  ! and atom_pairs_recv from the tasks lists
1076  CALL fb_com_tasks_build_atom_pairs(com_tasks=com_tasks_send, &
1077  atom_pairs=atom_pairs_send, &
1078  natoms_encode=nblkrows_total, &
1079  send_or_recv="send")
1080  CALL fb_com_tasks_build_atom_pairs(com_tasks=com_tasks_recv, &
1081  atom_pairs=atom_pairs_recv, &
1082  natoms_encode=nblkrows_total, &
1083  send_or_recv="recv")
1084 
1085  ! cleanup
1086  CALL fb_com_tasks_release(com_tasks_recv)
1087  CALL fb_com_tasks_release(com_tasks_send)
1088 
1089  CALL timestop(handle)
1090 
1091  END SUBROUTINE fb_fltrmat_generate_com_pairs
1092 
1093 ! **************************************************************************************************
1094 !> \brief generate the list of blocks (atom pairs) to be sent and received
1095 !> in order to construct the filter matrix for each atomic halo.
1096 !> This version is for use with fb_fltrmat_build_2, where MPI
1097 !> communications are done collectively.
1098 !> \param filter_mat : DBCSR formatted filter matrix
1099 !> \param atomic_halos : set of all local atomic halos contributing to the
1100 !> filter matrix
1101 !> \param para_env : cp2k parallel environment
1102 !> \param atom_pairs_send : list of blocks to be sent
1103 !> \param atom_pairs_recv : list of blocks to be received
1104 !> \author Lianheng Tong (LT) lianheng.tong@kcl.ac.uk
1105 ! **************************************************************************************************
1106  SUBROUTINE fb_fltrmat_generate_com_pairs_2(filter_mat, &
1107  atomic_halos, &
1108  para_env, &
1109  atom_pairs_send, &
1110  atom_pairs_recv)
1111  TYPE(dbcsr_type), POINTER :: filter_mat
1112  TYPE(fb_atomic_halo_list_obj), INTENT(IN) :: atomic_halos
1113  TYPE(mp_para_env_type), POINTER :: para_env
1114  TYPE(fb_com_atom_pairs_obj), INTENT(INOUT) :: atom_pairs_send, atom_pairs_recv
1115 
1116  CHARACTER(LEN=*), PARAMETER :: routinen = 'fb_fltrmat_generate_com_pairs_2'
1117 
1118  INTEGER :: dest, handle, iatom_global, iatom_in_halo, iatom_stored, ihalo, itask, &
1119  jatom_global, jatom_stored, natoms_in_halo, nblkrows_total, nhalos, ntasks_send
1120  INTEGER(KIND=int_8), DIMENSION(:, :), POINTER :: tasks_send
1121  INTEGER, DIMENSION(:), POINTER :: halo_atoms
1122  LOGICAL :: transpose
1123  TYPE(fb_atomic_halo_obj), DIMENSION(:), POINTER :: halos
1124  TYPE(fb_com_tasks_obj) :: com_tasks_recv, com_tasks_send
1125 
1126  CALL timeset(routinen, handle)
1127 
1128  NULLIFY (tasks_send)
1129  CALL fb_com_tasks_nullify(com_tasks_send)
1130  CALL fb_com_tasks_nullify(com_tasks_recv)
1131 
1132  ! initialise atom_pairs_send and atom_pairs_recv
1133  IF (fb_com_atom_pairs_has_data(atom_pairs_send)) THEN
1134  CALL fb_com_atom_pairs_init(atom_pairs_send)
1135  ELSE
1136  CALL fb_com_atom_pairs_create(atom_pairs_send)
1137  END IF
1138  IF (fb_com_atom_pairs_has_data(atom_pairs_recv)) THEN
1139  CALL fb_com_atom_pairs_init(atom_pairs_recv)
1140  ELSE
1141  CALL fb_com_atom_pairs_create(atom_pairs_recv)
1142  END IF
1143 
1144  ! The col atom block index for each filter matrix block are the
1145  ! owner atom of each halo. The row atom block index for each
1146  ! filter matrix block corresponding to each col are the halo atoms
1147  ! of the corresponding halos. Filter matrix is non-symmetric: it
1148  ! is non-square, because the blocks themselves are non-square
1149 
1150  CALL fb_atomic_halo_list_get(atomic_halos=atomic_halos, &
1151  nhalos=nhalos, &
1152  halos=halos)
1153 
1154  ! estimate the maximum number of blocks (i.e. atom paris) to send
1155  ntasks_send = 0
1156  DO ihalo = 1, nhalos
1157  CALL fb_atomic_halo_get(atomic_halo=halos(ihalo), &
1158  natoms=natoms_in_halo)
1159  ntasks_send = ntasks_send + natoms_in_halo
1160  END DO ! ihalo
1161 
1162  ! allocate send tasks
1163  ALLOCATE (tasks_send(task_n_records, ntasks_send))
1164 
1165  ! Get the total number of atoms. This can be obtained from the
1166  ! total number of block rows in the DBCSR filter matrix. We
1167  ! assumes that before calling this subroutine, the filter_mat has
1168  ! already been created and initialised: i.e. using
1169  ! dbcsr_create_new. Even if the matrix is at the moment empty,
1170  ! the attribute nblkrows_total is already assigned from the dbcsr
1171  ! distribution data
1172  CALL dbcsr_get_info(filter_mat, &
1173  nblkrows_total=nblkrows_total)
1174 
1175  ! source is always the local processor
1176  associate(src => para_env%mepos)
1177  ! construct send tasks
1178  itask = 1
1179  DO ihalo = 1, nhalos
1180  CALL fb_atomic_halo_get(atomic_halo=halos(ihalo), &
1181  owner_atom=jatom_global, &
1182  natoms=natoms_in_halo, &
1183  halo_atoms=halo_atoms)
1184  DO iatom_in_halo = 1, natoms_in_halo
1185  iatom_global = halo_atoms(iatom_in_halo)
1186  iatom_stored = iatom_global
1187  jatom_stored = jatom_global
1188  transpose = .false.
1189  ! find where the constructed block of filter matrix belongs to
1190  CALL dbcsr_get_stored_coordinates(filter_mat, &
1191  iatom_stored, &
1192  jatom_stored, &
1193  processor=dest)
1194  ! create the send tasks
1195  tasks_send(task_dest, itask) = dest
1196  tasks_send(task_src, itask) = src
1197  CALL fb_com_tasks_encode_pair(tasks_send(task_pair, itask), &
1198  iatom_global, jatom_global, &
1199  nblkrows_total)
1200  ! calculation of cost not implemented at the moment
1201  tasks_send(task_cost, itask) = 0
1202  itask = itask + 1
1203  END DO ! iatom_in_halo
1204  END DO ! ihalo
1205  END associate
1206 
1207  ! get the actual number of tasks
1208  ntasks_send = itask - 1
1209 
1210  CALL fb_com_tasks_create(com_tasks_send)
1211  CALL fb_com_tasks_set(com_tasks=com_tasks_send, &
1212  task_dim=task_n_records, &
1213  ntasks=ntasks_send, &
1214  nencode=nblkrows_total, &
1215  tasks=tasks_send)
1216 
1217  ! generate the recv task list (tasks_recv) from the send task list
1218  CALL fb_com_tasks_create(com_tasks_recv)
1219  CALL fb_com_tasks_transpose_dest_src(com_tasks_recv, "<", com_tasks_send, &
1220  para_env)
1221 
1222  ! task lists are now complete, now construct the atom_pairs_send
1223  ! and atom_pairs_recv from the tasks lists
1224  CALL fb_com_tasks_build_atom_pairs(com_tasks=com_tasks_send, &
1225  atom_pairs=atom_pairs_send, &
1226  natoms_encode=nblkrows_total, &
1227  send_or_recv="send")
1228  CALL fb_com_tasks_build_atom_pairs(com_tasks=com_tasks_recv, &
1229  atom_pairs=atom_pairs_recv, &
1230  natoms_encode=nblkrows_total, &
1231  send_or_recv="recv")
1232 
1233  ! cleanup
1234  CALL fb_com_tasks_release(com_tasks_recv)
1235  CALL fb_com_tasks_release(com_tasks_send)
1236 
1237  CALL timestop(handle)
1238 
1239  END SUBROUTINE fb_fltrmat_generate_com_pairs_2
1240 
1241 ! **************************************************************************************************
1242 !> \brief Build the atomic filter matrix for each atomic halo
1243 !> \param atomic_H : atomic KS matrix
1244 !> \param atomic_S : atomic overlap matrix
1245 !> \param fermi_level : fermi level used to construct the Fermi-Dirac
1246 !> filter function
1247 !> \param filter_temp : temperature used to construct the Fermi-Dirac
1248 !> filter function
1249 !> \param atomic_filter_mat : the atomic filter matrix
1250 !> \param tolerance : anything smaller than tolerance is treated as zero
1251 !> \author Lianheng Tong (LT) lianheng.tong@kcl.ac.uk
1252 ! **************************************************************************************************
1253  SUBROUTINE fb_fltrmat_build_atomic_fltrmat(atomic_H, &
1254  atomic_S, &
1255  fermi_level, &
1256  filter_temp, &
1257  atomic_filter_mat, &
1258  tolerance)
1259  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: atomic_h, atomic_s
1260  REAL(kind=dp), INTENT(IN) :: fermi_level, filter_temp
1261  REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: atomic_filter_mat
1262  REAL(kind=dp), INTENT(IN) :: tolerance
1263 
1264  CHARACTER(LEN=*), PARAMETER :: routinen = 'fb_fltrmat_build_atomic_fltrmat'
1265 
1266  INTEGER :: handle, handle_dgemm, handle_dsygv, ii, &
1267  info, jj, mat_dim, work_array_size
1268  LOGICAL :: check_ok
1269  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues, filter_function, work
1270  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: atomic_s_copy, eigenvectors, &
1271  filtered_eigenvectors
1272 
1273  CALL timeset(routinen, handle)
1274 
1275  ! This subroutine assumes atomic_filter_mat is not zero size, in
1276  ! other words, it really has to be constructed, instead of just
1277  ! being a dummy
1278 
1279  check_ok = SIZE(atomic_filter_mat, 1) > 0 .AND. &
1280  SIZE(atomic_filter_mat, 2) > 0
1281  cpassert(check_ok)
1282 
1283  ! initialise
1284  atomic_filter_mat = 0.0_dp
1285  mat_dim = SIZE(atomic_h, 1)
1286 
1287  ! diagonalise using LAPACK
1288  ALLOCATE (eigenvalues(mat_dim))
1289  ! get optimal work array size
1290  ALLOCATE (work(1))
1291  ! dsygv will overwrite part of atomic_H and atomic_S, thus need to copy them
1292  ALLOCATE (atomic_s_copy(SIZE(atomic_s, 1), SIZE(atomic_s, 2)))
1293  atomic_s_copy(:, :) = atomic_s(:, :)
1294  ALLOCATE (eigenvectors(SIZE(atomic_h, 1), SIZE(atomic_h, 2)))
1295  eigenvectors(:, :) = atomic_h(:, :)
1296 
1297  CALL timeset("fb_atomic_filter_dsygv", handle_dsygv)
1298 
1299  info = 0
1300  CALL dsygv(1, &
1301  'V', &
1302  'U', &
1303  mat_dim, &
1304  eigenvectors, &
1305  mat_dim, &
1306  atomic_s_copy, &
1307  mat_dim, &
1308  eigenvalues, &
1309  work, &
1310  -1, &
1311  info)
1312  work_array_size = nint(work(1))
1313  ! now allocate work array
1314  DEALLOCATE (work)
1315  ALLOCATE (work(work_array_size))
1316  work = 0.0_dp
1317  ! do calculation
1318  atomic_s_copy(:, :) = atomic_s(:, :)
1319  eigenvectors(:, :) = atomic_h(:, :)
1320  info = 0
1321  CALL dsygv(1, &
1322  'V', &
1323  'U', &
1324  mat_dim, &
1325  eigenvectors, &
1326  mat_dim, &
1327  atomic_s_copy, &
1328  mat_dim, &
1329  eigenvalues, &
1330  work, &
1331  work_array_size, &
1332  info)
1333  ! check if diagonalisation is successful
1334  IF (info .NE. 0) THEN
1335  WRITE (*, *) "DSYGV ERROR MESSAGE: ", info
1336  cpabort("DSYGV failed")
1337  END IF
1338 
1339  CALL timestop(handle_dsygv)
1340 
1341  DEALLOCATE (work)
1342  DEALLOCATE (atomic_s_copy)
1343 
1344  ! first get the filter function
1345  ALLOCATE (filter_function(mat_dim))
1346  filter_function = 0.0_dp
1347  CALL fb_fltrmat_fermi_dirac_mu(filter_function, &
1348  eigenvalues, &
1349  filter_temp, &
1350  fermi_level)
1351  DEALLOCATE (eigenvalues)
1352 
1353  ! atomic_H has the eigenvectors, construct the version of it
1354  ! filtered through the filter function
1355  ALLOCATE (filtered_eigenvectors(mat_dim, mat_dim))
1356  DO jj = 1, mat_dim
1357  DO ii = 1, mat_dim
1358  filtered_eigenvectors(ii, jj) = &
1359  filter_function(jj)*eigenvectors(ii, jj)
1360  END DO ! ii
1361  END DO ! jj
1362 
1363  DEALLOCATE (filter_function)
1364 
1365  CALL timeset("fb_atomic_filter_dgemm", handle_dgemm)
1366 
1367  ! construct atomic filter matrix
1368  CALL dgemm("N", &
1369  "T", &
1370  mat_dim, &
1371  mat_dim, &
1372  mat_dim, &
1373  1.0_dp, &
1374  filtered_eigenvectors, &
1375  mat_dim, &
1376  eigenvectors, &
1377  mat_dim, &
1378  0.0_dp, &
1379  atomic_filter_mat, &
1380  mat_dim)
1381 
1382  CALL timestop(handle_dgemm)
1383 
1384  ! remove small negative terms due to numerical error, the filter
1385  ! matrix must not be negative definite
1386  DO jj = 1, SIZE(atomic_filter_mat, 2)
1387  DO ii = 1, SIZE(atomic_filter_mat, 1)
1388  IF (abs(atomic_filter_mat(ii, jj)) < tolerance) THEN
1389  atomic_filter_mat(ii, jj) = 0.0_dp
1390  END IF
1391  END DO
1392  END DO
1393 
1394  DEALLOCATE (filtered_eigenvectors)
1395  DEALLOCATE (eigenvectors)
1396 
1397  CALL timestop(handle)
1398 
1399  END SUBROUTINE fb_fltrmat_build_atomic_fltrmat
1400 
1401 ! **************************************************************************************************
1402 !> \brief get values of Fermi-Dirac distribution based on a given fermi
1403 !> level at a given set of energy eigenvalues
1404 !> \param f : the Fermi-Dirac distribution function values
1405 !> \param eigenvals : set of energy eigenvalues
1406 !> \param T : temperature
1407 !> \param mu : the fermi level
1408 !> \author Lianheng Tong (LT) lianheng.tong@kcl.ac.uk
1409 ! **************************************************************************************************
1410  SUBROUTINE fb_fltrmat_fermi_dirac_mu(f, eigenvals, T, mu)
1411  REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: f
1412  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: eigenvals
1413  REAL(kind=dp), INTENT(IN) :: t, mu
1414 
1415  REAL(kind=dp) :: kts, ne
1416 
1417 ! we want fermi function max at 1, so maxocc = 1 here
1418 
1419  CALL fermi(f, ne, kts, eigenvals, mu, t, 1.0_dp)
1420  END SUBROUTINE fb_fltrmat_fermi_dirac_mu
1421 
1422 ! **************************************************************************************************
1423 !> \brief get values of Fermi-Dirac distribution based on a given electron
1424 !> number at a given set of energy eigenvales
1425 !> \param f : the Fermi-Dirac distribution function values
1426 !> \param eigenvals : set of energy eigenvalues
1427 !> \param T : temperature
1428 !> \param ne : number of electrons
1429 !> \param maxocc : maximum occupancy per orbital
1430 !> \author Lianheng Tong (LT) lianheng.tong@kcl.ac.uk
1431 ! **************************************************************************************************
1432  SUBROUTINE fb_fltrmat_fermi_dirac_ne(f, eigenvals, T, ne, maxocc)
1433  REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: f
1434  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: eigenvals
1435  REAL(kind=dp), INTENT(IN) :: t, ne, maxocc
1436 
1437  REAL(kind=dp) :: kts, mu
1438 
1439 ! mu is the calculated fermi level
1440 ! kTS is the calculated entropic contribution to the energy i.e. -TS
1441 ! kTS = kT*[f ln f + (1-f) ln (1-f)]
1442 
1443  CALL fermifixed(f, mu, kts, eigenvals, ne, t, maxocc)
1444  END SUBROUTINE fb_fltrmat_fermi_dirac_ne
1445 
1446 END MODULE qs_fb_filter_matrix_methods
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
deal with the Fermi distribution, compute it, fix mu, get derivs
Definition: fermi_utils.F:13
subroutine, public fermi(f, N, kTS, e, mu, T, maxocc, estate, festate)
returns occupations according to Fermi-Dirac statistics for a given set of energies and fermi level....
Definition: fermi_utils.F:48
subroutine, public fermifixed(f, mu, kTS, e, N, T, maxocc, estate, festate)
returns occupations according to Fermi-Dirac statistics for a given set of energies and number of ele...
Definition: fermi_utils.F:202
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public int_8
Definition: kinds.F:54
integer, parameter, public dp
Definition: kinds.F:34
integer, parameter, public default_string_length
Definition: kinds.F:57
Interface to the message passing library MPI.
Define the data structure for the particle information.
subroutine, public fb_atomic_halo_list_get(atomic_halos, nhalos, max_nhalos, halos)
Gets attributes from an fb_atomic_halo_list object, one should only access the data content in a fb_a...
subroutine, public fb_atomic_halo_get(atomic_halo, owner_atom, owner_id_in_halo, natoms, nelectrons, halo_atoms, sorted, cost)
Gets attributes from a fb_atomic_halo object, one should only access the data content in a fb_atomic_...
subroutine, public fb_atomic_halo_set(atomic_halo, owner_atom, owner_id_in_halo, natoms, nelectrons, halo_atoms, sorted, cost)
Sets attributes in a fb_atomic_halo object, one should only set the data content in a fb_atomic_halo ...
subroutine, public fb_atomic_halo_create(atomic_halo)
Creates and initialises an empty fb_atomic_halo object.
subroutine, public fb_atomic_halo_nullify(atomic_halo)
Nullifies a fb_atomic_halo object, note that it does not release the original object....
subroutine, public fb_atomic_halo_release(atomic_halo)
Releases an fb_atomic_halo object.
subroutine, public fb_atmatrix_calc_size(dbcsr_mat, atomic_halo, nrows, ncols, blk_row_start, blk_col_start)
Calculates the atomic matrix size from a given DBCSR matrix and atomic halo. It also calculates the f...
subroutine, public fb_atmatrix_generate_com_pairs_2(dbcsr_mat, atomic_halos, para_env, atom_pairs_send, atom_pairs_recv)
generate list of blocks (atom pairs) of a DBCSR matrix to be sent and received in order to construct ...
subroutine, public fb_atmatrix_construct_2(matrix_storage, atomic_halo, atomic_matrix, blk_row_start, blk_col_start)
Constructs atomic matrix for filter basis method from a given DBCSR matrix and a set of atomic send a...
subroutine, public fb_atmatrix_construct(dbcsr_mat, atomic_halo, para_env, atomic_matrix, blk_row_start, blk_col_start)
Constructs atomic matrix for filter basis method from a given DBCSR matrix and a set of atomic send a...
subroutine, public fb_com_atom_pairs_init(atom_pairs)
Initialises an fb_com_atom_pairs object, and makes it empty.
integer, parameter, public task_pair
subroutine, public fb_com_atom_pairs_decode(ind, pe, iatom, jatom, natoms)
Decodes a single integer into the (rank, iatom, jatom) index of a communication task to send/receive ...
integer, parameter, public task_src
subroutine, public fb_com_tasks_nullify(com_tasks)
Nullifies a fb_com_tasks object, note that it does not release the original object....
subroutine, public fb_com_tasks_encode_pair(ind, iatom, jatom, natoms)
Encodes (iatom, jatom) pair index of a block into a single integer.
subroutine, public fb_com_atom_pairs_distribute_blks(matrix_storage, atom_pairs_send, atom_pairs_recv, para_env, dbcsr_mat)
Given send and recv fb_com_atom_pair object, distribute the matrix blocks stored in a fb_matrix_data ...
subroutine, public fb_com_tasks_build_atom_pairs(com_tasks, atom_pairs, natoms_encode, send_or_recv)
Generate send or receive atom_pair lists from a com_tasks object. atom_pair list is used as a condens...
subroutine, public fb_com_atom_pairs_create(atom_pairs)
Creates and initialises an empty fb_com_atom_pairs object.
subroutine, public fb_com_tasks_transpose_dest_src(tasks_dest_is_me, direction, tasks_src_is_me, para_env)
Start from a local set of tasks that has desc/src process equal to the local MPI rank,...
logical function, public fb_com_atom_pairs_has_data(atom_pairs)
Checks if a fb_com_atom_pairs object is associated with an actual data content or not.
subroutine, public fb_com_atom_pairs_get(atom_pairs, npairs, natoms_encode, pairs)
Gets attributes from a fb_com_atom_pairs object, one should only access the data content in a fb_com_...
subroutine, public fb_com_atom_pairs_nullify(atom_pairs)
Nullifies a fb_com_atom_pairs object, note that it does not release the original object....
subroutine, public fb_com_atom_pairs_release(atom_pairs)
Releases an fb_com_atom_pairs object.
subroutine, public fb_com_tasks_create(com_tasks)
Creates and initialises an empty fb_com_tasks object.
integer, parameter, public task_cost
subroutine, public fb_com_tasks_release(com_tasks)
Releases an fb_com_tasks object.
subroutine, public fb_com_atom_pairs_gather_blks(dbcsr_mat, atom_pairs_send, atom_pairs_recv, para_env, matrix_storage)
Given send and recv fb_com_atom_pair object, gather all the relevant DBCSR matrix blocks together,...
subroutine, public fb_com_tasks_set(com_tasks, task_dim, ntasks, nencode, tasks)
Sets attributes in a fb_com_tasks object, one should only access the data content in a fb_com_tasks o...
integer, parameter, public task_dest
integer, parameter, public task_n_records
subroutine, public fb_com_atom_pairs_calc_buffer_sizes(atom_pairs, nprocs, row_blk_sizes, col_blk_sizes, sendrecv_sizes, sendrecv_disps, sendrecv_pair_counts, sendrecv_pair_disps)
Calculate the MPI send or recv buffer sizes according to the communication pairs (atom_pairs) and DBC...
subroutine, public fb_fltrmat_build_2(H_mat, S_mat, atomic_halos, trial_fns, para_env, particle_set, fermi_level, filter_temp, name, filter_mat, tolerance)
Build the filter matrix, with MPI communications grouped together. More effcient on communication,...
subroutine, public fb_fltrmat_build(H_mat, S_mat, atomic_halos, trial_fns, para_env, particle_set, fermi_level, filter_temp, name, filter_mat, tolerance)
Build the filter matrix, with MPI communications happening at each step. Less efficient on communicat...
pure logical function, public fb_matrix_data_has_data(matrix_data)
check if the object has data associated to it
subroutine, public fb_matrix_data_add(matrix_data, row, col, blk)
Add a matrix block to a fb_matrix_data object.
subroutine, public fb_matrix_data_create(matrix_data, nmax, nencode)
Creates and initialises an empty fb_matrix_data object of a given size.
pure subroutine, public fb_matrix_data_nullify(matrix_data)
Nullifies a fb_matrix_data object.
subroutine, public fb_matrix_data_release(matrix_data)
releases given object
subroutine, public fb_trial_fns_get(trial_fns, nfunctions, functions)
get values of the attributes of a fb_trial_fns object
Utilities for string manipulations.
subroutine, public compress(string, full)
Eliminate multiple space characters in a string. If full is .TRUE., then all spaces are eliminated.
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.