LCOV - code coverage report
Current view: top level - src - qs_fb_filter_matrix_methods.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:71c3ab0) Lines: 97.3 % 403 392
Test Date: 2026-07-25 06:35:44 Functions: 88.9 % 9 8

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : MODULE qs_fb_filter_matrix_methods
       9              : 
      10              :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      11              :                                               get_atomic_kind
      12              :    USE cp_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 input_constants,                 ONLY: smear_fermi_dirac
      21              :    USE kinds,                           ONLY: default_string_length,&
      22              :                                               dp,&
      23              :                                               int_8
      24              :    USE message_passing,                 ONLY: mp_para_env_type
      25              :    USE particle_types,                  ONLY: particle_type
      26              :    USE qs_fb_atomic_halo_types,         ONLY: fb_atomic_halo_create,&
      27              :                                               fb_atomic_halo_get,&
      28              :                                               fb_atomic_halo_list_get,&
      29              :                                               fb_atomic_halo_list_obj,&
      30              :                                               fb_atomic_halo_nullify,&
      31              :                                               fb_atomic_halo_obj,&
      32              :                                               fb_atomic_halo_release,&
      33              :                                               fb_atomic_halo_set
      34              :    USE qs_fb_atomic_matrix_methods,     ONLY: fb_atmatrix_calc_size,&
      35              :                                               fb_atmatrix_construct,&
      36              :                                               fb_atmatrix_construct_2,&
      37              :                                               fb_atmatrix_generate_com_pairs_2
      38              :    USE qs_fb_com_tasks_types,           ONLY: &
      39              :         TASK_COST, TASK_DEST, TASK_N_RECORDS, TASK_PAIR, TASK_SRC, &
      40              :         fb_com_atom_pairs_calc_buffer_sizes, fb_com_atom_pairs_create, fb_com_atom_pairs_decode, &
      41              :         fb_com_atom_pairs_distribute_blks, fb_com_atom_pairs_gather_blks, fb_com_atom_pairs_get, &
      42              :         fb_com_atom_pairs_has_data, fb_com_atom_pairs_init, fb_com_atom_pairs_nullify, &
      43              :         fb_com_atom_pairs_obj, fb_com_atom_pairs_release, fb_com_tasks_build_atom_pairs, &
      44              :         fb_com_tasks_create, fb_com_tasks_encode_pair, fb_com_tasks_nullify, fb_com_tasks_obj, &
      45              :         fb_com_tasks_release, fb_com_tasks_set, fb_com_tasks_transpose_dest_src
      46              :    USE qs_fb_matrix_data_types,         ONLY: fb_matrix_data_add,&
      47              :                                               fb_matrix_data_create,&
      48              :                                               fb_matrix_data_has_data,&
      49              :                                               fb_matrix_data_nullify,&
      50              :                                               fb_matrix_data_obj,&
      51              :                                               fb_matrix_data_release
      52              :    USE qs_fb_trial_fns_types,           ONLY: fb_trial_fns_get,&
      53              :                                               fb_trial_fns_obj
      54              :    USE smearing_utils,                  ONLY: SmearFixed,&
      55              :                                               SmearOcc
      56              :    USE string_utilities,                ONLY: compress,&
      57              :                                               uppercase
      58              : #include "./base/base_uses.f90"
      59              : 
      60              :    IMPLICIT NONE
      61              : 
      62              :    PRIVATE
      63              : 
      64              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_fb_filter_matrix_methods'
      65              : 
      66              :    PUBLIC :: fb_fltrmat_build, &
      67              :              fb_fltrmat_build_2
      68              : 
      69              : CONTAINS
      70              : 
      71              : ! **************************************************************************************************
      72              : !> \brief Build the filter matrix, with MPI communications happening at each
      73              : !>        step. Less efficient on communication, but more efficient on
      74              : !>        memory usage (compared to fb_fltrmat_build_2)
      75              : !> \param H_mat : DBCSR system KS matrix
      76              : !> \param S_mat : DBCSR system overlap matrix
      77              : !> \param atomic_halos : list of all local atomic halos, each halo gives
      78              : !>                       one atomic matrix and contributes to one blk
      79              : !>                       col to the filter matrix
      80              : !> \param trial_fns : the trial functions to be used to shrink the
      81              : !>                     size of the new "filtered" basis
      82              : !> \param para_env : cp2k parallel environment
      83              : !> \param particle_set : set of all particles in the system
      84              : !> \param fermi_level : the fermi level used for defining the filter
      85              : !>                      function, which is a Fermi-Dirac distribution
      86              : !>                      function
      87              : !> \param filter_temp : the filter temperature used for defining the
      88              : !>                      filter function
      89              : !> \param name        : name given to the filter matrix
      90              : !> \param filter_mat  : DBCSR format filter matrix
      91              : !> \param tolerance   : anything less than tolerance is treated as zero
      92              : !> \author Lianheng Tong (LT) lianheng.tong@kcl.ac.uk
      93              : ! **************************************************************************************************
      94           64 :    SUBROUTINE fb_fltrmat_build(H_mat, &
      95              :                                S_mat, &
      96              :                                atomic_halos, &
      97              :                                trial_fns, &
      98              :                                para_env, &
      99              :                                particle_set, &
     100              :                                fermi_level, &
     101              :                                filter_temp, &
     102              :                                name, &
     103              :                                filter_mat, &
     104              :                                tolerance)
     105              :       TYPE(dbcsr_type), POINTER                          :: H_mat, S_mat
     106              :       TYPE(fb_atomic_halo_list_obj), INTENT(IN)          :: atomic_halos
     107              :       TYPE(fb_trial_fns_obj), INTENT(IN)                 :: trial_fns
     108              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     109              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     110              :       REAL(KIND=dp), INTENT(IN)                          :: fermi_level, filter_temp
     111              :       CHARACTER(LEN=*), INTENT(IN)                       :: name
     112              :       TYPE(dbcsr_type), POINTER                          :: filter_mat
     113              :       REAL(KIND=dp), INTENT(IN)                          :: tolerance
     114              : 
     115              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'fb_fltrmat_build'
     116              : 
     117              :       CHARACTER(LEN=32)                                  :: symmetry_string
     118              :       CHARACTER(LEN=default_string_length)               :: name_string
     119              :       INTEGER                                            :: handle, iblkcol, ihalo, ikind, &
     120              :                                                             max_nhalos, nblkcols_total, nhalos
     121           64 :       INTEGER, DIMENSION(:), POINTER                     :: col_blk_size, dummy_halo_atoms, ntfns, &
     122           64 :                                                             row_blk_size
     123              :       LOGICAL                                            :: send_data_only
     124              :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
     125              :       TYPE(dbcsr_distribution_type)                      :: dbcsr_dist
     126              :       TYPE(fb_atomic_halo_obj)                           :: dummy_atomic_halo
     127           64 :       TYPE(fb_atomic_halo_obj), DIMENSION(:), POINTER    :: halos
     128              : 
     129           64 :       CALL timeset(routineN, handle)
     130              : 
     131           64 :       NULLIFY (halos, atomic_kind, ntfns, dummy_halo_atoms, row_blk_size, col_blk_size)
     132           64 :       CALL fb_atomic_halo_nullify(dummy_atomic_halo)
     133              : 
     134              :       ! filter_mat must be of a dissassociated status (i.e. brand new)
     135           64 :       CPASSERT(.NOT. ASSOCIATED(filter_mat))
     136              : 
     137              :       ! get trial function information
     138              :       CALL fb_trial_fns_get(trial_fns=trial_fns, &
     139           64 :                             nfunctions=ntfns)
     140              : 
     141              :       ! calculate the row_blk_size and col_blk_size arrays for
     142              :       ! constructing the filter matrix in DBCSR format
     143              :       ! row_blk_size for the filter matrix is the same as H or S
     144              :       CALL dbcsr_get_info(H_mat, &
     145              :                           nblkcols_total=nblkcols_total, &
     146              :                           row_blk_size=row_blk_size, &
     147           64 :                           distribution=dbcsr_dist)
     148          192 :       ALLOCATE (col_blk_size(nblkcols_total))
     149          576 :       col_blk_size = 0
     150          576 :       DO iblkcol = 1, nblkcols_total
     151          512 :          atomic_kind => particle_set(iblkcol)%atomic_kind
     152              :          CALL get_atomic_kind(atomic_kind=atomic_kind, &
     153          512 :                               kind_number=ikind)
     154          576 :          col_blk_size(iblkcol) = ntfns(ikind)
     155              :       END DO
     156              :       ! DO NOT deallocate cbs if gift=.TRUE. as col_blk_sizes will only point to cbs
     157           64 :       name_string = name
     158           64 :       CALL compress(name_string)
     159           64 :       CALL uppercase(name_string)
     160              :       ! the filter matrix is non-square and is always non-symmetric
     161           64 :       symmetry_string = dbcsr_type_no_symmetry
     162              :       ! create empty filter matrix
     163           64 :       ALLOCATE (filter_mat)
     164              :       CALL dbcsr_create(matrix=filter_mat, &
     165              :                         name=name_string, &
     166              :                         dist=dbcsr_dist, &
     167              :                         matrix_type=symmetry_string, &
     168              :                         row_blk_size=row_blk_size, &
     169           64 :                         col_blk_size=col_blk_size)
     170           64 :       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           64 :                                    halos=halos)
     176              : 
     177              :       ! create dummy empty atomic halo
     178           64 :       CALL fb_atomic_halo_create(dummy_atomic_halo)
     179           64 :       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           64 :                               sorted=.TRUE.)
     187              : 
     188           64 :       send_data_only = .FALSE.
     189              : 
     190          320 :       DO ihalo = 1, max_nhalos
     191          256 :          IF (ihalo > nhalos) THEN
     192              :             send_data_only = .TRUE.
     193              :          END IF
     194              :          ! construct the filter matrix block by block
     195          320 :          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            0 :                                        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          256 :                                        tolerance)
     217              :          END IF ! send_data_only
     218              :       END DO
     219              : 
     220              :       ! finalise the filter matrix
     221           64 :       CALL dbcsr_finalize(filter_mat)
     222              : 
     223              :       ! cleanup
     224           64 :       CALL fb_atomic_halo_release(dummy_atomic_halo)
     225              : 
     226           64 :       CALL timestop(handle)
     227              : 
     228          192 :    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           16 :    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           16 :       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           16 :       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           16 :       CALL timeset(routineN, handle)
     292              : 
     293           16 :       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           16 :       check_ok = .NOT. ASSOCIATED(filter_mat)
     297           16 :       CPASSERT(check_ok)
     298              : 
     299              :       ! get total number of atoms
     300           16 :       natoms_global = SIZE(particle_set)
     301              : 
     302              :       ! get trial function information
     303              :       CALL fb_trial_fns_get(trial_fns=trial_fns, &
     304           16 :                             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           16 :                           distribution=dbcsr_dist)
     313           48 :       ALLOCATE (col_blk_size(nblkcols_total))
     314          144 :       col_blk_size = 0
     315          144 :       DO iblkcol = 1, nblkcols_total
     316          128 :          atomic_kind => particle_set(iblkcol)%atomic_kind
     317              :          CALL get_atomic_kind(atomic_kind=atomic_kind, &
     318          128 :                               kind_number=ikind)
     319          144 :          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           16 :       name_string = name
     323           16 :       CALL compress(name_string)
     324           16 :       CALL uppercase(name_string)
     325              :       ! create empty filter matrix (it is always non-symmetric as it is non-square)
     326           16 :       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           16 :                         col_blk_size=col_blk_size)
     333           16 :       DEALLOCATE (col_blk_size)
     334              : 
     335              :       ! get all the blocks required for constructing atomic matrics, and
     336              :       ! store it in a fb_matrix_data object
     337           16 :       CALL fb_matrix_data_nullify(H_mat_data)
     338           16 :       CALL fb_matrix_data_nullify(S_mat_data)
     339           16 :       CALL fb_com_atom_pairs_nullify(atmatrix_blks_send)
     340           16 :       CALL fb_com_atom_pairs_nullify(atmatrix_blks_recv)
     341           16 :       CALL fb_com_atom_pairs_create(atmatrix_blks_send)
     342           16 :       CALL fb_com_atom_pairs_create(atmatrix_blks_recv)
     343              :       ! H matrix
     344              :       CALL fb_atmatrix_generate_com_pairs_2(H_mat, &
     345              :                                             atomic_halos, &
     346              :                                             para_env, &
     347              :                                             atmatrix_blks_send, &
     348           16 :                                             atmatrix_blks_recv)
     349              :       CALL fb_com_atom_pairs_get(atom_pairs=atmatrix_blks_recv, &
     350           16 :                                  npairs=nblks_recv)
     351              :       CALL fb_matrix_data_create(H_mat_data, &
     352              :                                  nblks_recv, &
     353           16 :                                  natoms_global)
     354              :       CALL fb_com_atom_pairs_gather_blks(H_mat, &
     355              :                                          atmatrix_blks_send, &
     356              :                                          atmatrix_blks_recv, &
     357              :                                          para_env, &
     358           16 :                                          H_mat_data)
     359              :       ! S matrix
     360              :       CALL fb_atmatrix_generate_com_pairs_2(S_mat, &
     361              :                                             atomic_halos, &
     362              :                                             para_env, &
     363              :                                             atmatrix_blks_send, &
     364           16 :                                             atmatrix_blks_recv)
     365              :       CALL fb_com_atom_pairs_get(atom_pairs=atmatrix_blks_recv, &
     366           16 :                                  npairs=nblks_recv)
     367              :       CALL fb_matrix_data_create(S_mat_data, &
     368              :                                  nblks_recv, &
     369           16 :                                  natoms_global)
     370              :       CALL fb_com_atom_pairs_gather_blks(S_mat, &
     371              :                                          atmatrix_blks_send, &
     372              :                                          atmatrix_blks_recv, &
     373              :                                          para_env, &
     374           16 :                                          S_mat_data)
     375              :       ! cleanup
     376           16 :       CALL fb_com_atom_pairs_release(atmatrix_blks_send)
     377           16 :       CALL fb_com_atom_pairs_release(atmatrix_blks_recv)
     378              : 
     379              :       ! make filter matrix blocks one by one and store in an
     380              :       ! matrix_data_obj
     381           16 :       CALL fb_matrix_data_nullify(filter_mat_data)
     382              :       CALL fb_atomic_halo_list_get(atomic_halos=atomic_halos, &
     383              :                                    nhalos=nhalos, &
     384           16 :                                    halos=halos)
     385           16 :       nmax = 0
     386           80 :       DO ihalo = 1, nhalos
     387              :          CALL fb_atomic_halo_get(atomic_halo=halos(ihalo), &
     388           64 :                                  natoms=natoms_in_halo)
     389           80 :          nmax = nmax + natoms_in_halo
     390              :       END DO
     391              :       CALL fb_matrix_data_create(filter_mat_data, &
     392              :                                  nmax, &
     393           16 :                                  natoms_global)
     394           80 :       DO ihalo = 1, nhalos
     395              :          CALL fb_fltrmat_add_blkcol_2(H_mat, &
     396              :                                       S_mat, &
     397              :                                       H_mat_data, &
     398              :                                       S_mat_data, &
     399              :                                       halos(ihalo), &
     400              :                                       trial_fns, &
     401              :                                       particle_set, &
     402              :                                       fermi_level, &
     403              :                                       filter_temp, &
     404              :                                       filter_mat_data, &
     405           80 :                                       tolerance)
     406              :       END DO
     407              :       ! clean up
     408           16 :       CALL fb_matrix_data_release(H_mat_data)
     409           16 :       CALL fb_matrix_data_release(S_mat_data)
     410              : 
     411              :       ! distribute the relevant blocks from the matrix_data_obj to DBCSR
     412              :       ! filter matrix
     413           16 :       CALL fb_com_atom_pairs_nullify(filter_mat_blks_send)
     414           16 :       CALL fb_com_atom_pairs_nullify(filter_mat_blks_recv)
     415           16 :       CALL fb_com_atom_pairs_create(filter_mat_blks_send)
     416           16 :       CALL fb_com_atom_pairs_create(filter_mat_blks_recv)
     417              :       CALL fb_fltrmat_generate_com_pairs_2(filter_mat, &
     418              :                                            atomic_halos, &
     419              :                                            para_env, &
     420              :                                            filter_mat_blks_send, &
     421           16 :                                            filter_mat_blks_recv)
     422              :       CALL fb_com_atom_pairs_distribute_blks(filter_mat_data, &
     423              :                                              filter_mat_blks_send, &
     424              :                                              filter_mat_blks_recv, &
     425              :                                              para_env, &
     426           16 :                                              filter_mat)
     427              :       ! cleanup
     428           16 :       CALL fb_com_atom_pairs_release(filter_mat_blks_send)
     429           16 :       CALL fb_com_atom_pairs_release(filter_mat_blks_recv)
     430           16 :       CALL fb_matrix_data_release(filter_mat_data)
     431              : 
     432              :       ! finalise matrix
     433           16 :       CALL dbcsr_finalize(filter_mat)
     434              : 
     435           16 :       CALL timestop(handle)
     436              : 
     437           96 :    END SUBROUTINE fb_fltrmat_build_2
     438              : 
     439              : ! **************************************************************************************************
     440              : !> \brief Add a computed blocks in one column to the filter matrix. This
     441              : !>        version is used by fb_fltrmat_build, for the case where MPI
     442              : !>        communications are done at each step
     443              : !>        It does not finalise the filter matrix
     444              : !> \param H_mat : DBCSR system KS matrix
     445              : !> \param S_mat : DBCSR system overlap matrix
     446              : !> \param atomic_halo :  the halo that contributes to the blk
     447              : !>                       col of the filter matrix
     448              : !> \param trial_fns ...
     449              : !> \param para_env : cp2k parallel environment
     450              : !> \param particle_set : set of all particles in the system
     451              : !> \param fermi_level : the fermi level used for defining the filter
     452              : !>                      function, which is a Fermi-Dirac distribution
     453              : !>                      function
     454              : !> \param filter_temp : the filter temperature used for defining the
     455              : !>                      filter function
     456              : !> \param filter_mat  : DBCSR format filter matrix
     457              : !> \param tolerance   : anything smaller than tolerance is treated as zero
     458              : !> \author Lianheng Tong (LT) lianheng.tong@kcl.ac.uk
     459              : ! **************************************************************************************************
     460          256 :    SUBROUTINE fb_fltrmat_add_blkcol(H_mat, &
     461              :                                     S_mat, &
     462              :                                     atomic_halo, &
     463              :                                     trial_fns, &
     464              :                                     para_env, &
     465              :                                     particle_set, &
     466              :                                     fermi_level, &
     467              :                                     filter_temp, &
     468              :                                     filter_mat, &
     469              :                                     tolerance)
     470              :       TYPE(dbcsr_type), POINTER                          :: H_mat, S_mat
     471              :       TYPE(fb_atomic_halo_obj), INTENT(IN)               :: atomic_halo
     472              :       TYPE(fb_trial_fns_obj), INTENT(IN)                 :: trial_fns
     473              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     474              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     475              :       REAL(KIND=dp), INTENT(IN)                          :: fermi_level, filter_temp
     476              :       TYPE(dbcsr_type), POINTER                          :: filter_mat
     477              :       REAL(KIND=dp), INTENT(IN)                          :: tolerance
     478              : 
     479              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'fb_fltrmat_add_blkcol'
     480              : 
     481              :       INTEGER :: handle, handle_mpi, iatom_global, iatom_in_halo, ind, ipair, ipe, itrial, &
     482              :          jatom_global, jatom_in_halo, jkind, natoms_global, natoms_in_halo, ncols_atmatrix, &
     483              :          ncols_blk, nrows_atmatrix, nrows_blk, numprocs, pe, recv_encode, send_encode
     484          256 :       INTEGER(KIND=int_8), DIMENSION(:), POINTER         :: pairs_recv, pairs_send
     485          256 :       INTEGER, ALLOCATABLE, DIMENSION(:) :: atomic_H_blk_col_start, atomic_H_blk_row_start, &
     486          256 :          atomic_S_blk_col_start, atomic_S_blk_row_start, col_block_size_data, ind_in_halo, &
     487          256 :          recv_disps, recv_pair_count, recv_pair_disps, recv_sizes, send_disps, send_pair_count, &
     488          256 :          send_pair_disps, send_sizes
     489          256 :       INTEGER, DIMENSION(:), POINTER                     :: halo_atoms, ntfns, row_block_size_data
     490          256 :       INTEGER, DIMENSION(:, :), POINTER                  :: tfns
     491          256 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), TARGET   :: recv_buf, send_buf
     492          256 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: atomic_filter_mat, atomic_H, atomic_S
     493          256 :       REAL(kind=dp), DIMENSION(:), POINTER               :: vector
     494              :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
     495              :       TYPE(fb_com_atom_pairs_obj)                        :: com_pairs_recv, com_pairs_send
     496              : 
     497          256 :       CALL timeset(routineN, handle)
     498              : 
     499          256 :       NULLIFY (atomic_kind, halo_atoms, ntfns, pairs_send, pairs_recv, &
     500          256 :                row_block_size_data, tfns)
     501          256 :       CALL fb_com_atom_pairs_nullify(com_pairs_send)
     502          256 :       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          256 :       CALL fb_com_atom_pairs_create(com_pairs_send)
     510          256 :       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          256 :                                          com_pairs_recv)
     516              :       CALL fb_com_atom_pairs_get(atom_pairs=com_pairs_send, &
     517              :                                  natoms_encode=send_encode, &
     518          256 :                                  pairs=pairs_send)
     519              :       CALL fb_com_atom_pairs_get(atom_pairs=com_pairs_recv, &
     520              :                                  natoms_encode=recv_encode, &
     521          256 :                                  pairs=pairs_recv)
     522              : 
     523              :       ! get para_env info
     524          256 :       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          256 :                             functions=tfns)
     531              : 
     532              :       ! obtain row and col block size data for filter matrix
     533          256 :       CALL dbcsr_get_info(H_mat, row_blk_size=row_block_size_data)
     534          256 :       natoms_global = SIZE(particle_set)
     535          768 :       ALLOCATE (col_block_size_data(natoms_global))
     536         2304 :       DO jatom_global = 1, natoms_global
     537         2048 :          atomic_kind => particle_set(jatom_global)%atomic_kind
     538         2048 :          CALL get_atomic_kind(atomic_kind=atomic_kind, kind_number=jkind)
     539         2304 :          col_block_size_data(jatom_global) = ntfns(jkind)
     540              :       END DO
     541              : 
     542              :       ! allocate temporary arrays for send
     543          768 :       ALLOCATE (send_sizes(numprocs))
     544          512 :       ALLOCATE (send_disps(numprocs))
     545          512 :       ALLOCATE (send_pair_count(numprocs))
     546          512 :       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          256 :                                                send_pair_disps)
     556              :       ! allocate send buffer
     557         1280 :       ALLOCATE (send_buf(SUM(send_sizes)))
     558              : 
     559              :       ! allocate temporary array for recv
     560          512 :       ALLOCATE (recv_sizes(numprocs))
     561          512 :       ALLOCATE (recv_disps(numprocs))
     562          512 :       ALLOCATE (recv_pair_count(numprocs))
     563          512 :       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          256 :                                                recv_pair_disps)
     573              :       ! allocate recv buffer
     574         1280 :       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          256 :                               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         1024 :                 atomic_H_blk_col_start(natoms_in_halo + 1))
     587              :       CALL fb_atmatrix_calc_size(H_mat, &
     588              :                                  atomic_halo, &
     589              :                                  nrows_atmatrix, &
     590              :                                  ncols_atmatrix, &
     591              :                                  atomic_H_blk_row_start, &
     592          256 :                                  atomic_H_blk_col_start)
     593              : 
     594         1024 :       ALLOCATE (atomic_H(nrows_atmatrix, ncols_atmatrix))
     595              :       CALL fb_atmatrix_construct(H_mat, &
     596              :                                  atomic_halo, &
     597              :                                  para_env, &
     598              :                                  atomic_H, &
     599              :                                  atomic_H_blk_row_start, &
     600          256 :                                  atomic_H_blk_col_start)
     601              : 
     602              :       ! construct atomic matrix for S for atomic_halo
     603              :       ALLOCATE (atomic_S_blk_row_start(natoms_in_halo + 1), &
     604          768 :                 atomic_S_blk_col_start(natoms_in_halo + 1))
     605              :       CALL fb_atmatrix_calc_size(S_mat, &
     606              :                                  atomic_halo, &
     607              :                                  nrows_atmatrix, &
     608              :                                  ncols_atmatrix, &
     609              :                                  atomic_S_blk_row_start, &
     610          256 :                                  atomic_S_blk_col_start)
     611         1024 :       ALLOCATE (atomic_S(nrows_atmatrix, ncols_atmatrix))
     612              :       CALL fb_atmatrix_construct(S_mat, &
     613              :                                  atomic_halo, &
     614              :                                  para_env, &
     615              :                                  atomic_S, &
     616              :                                  atomic_S_blk_row_start, &
     617          256 :                                  atomic_S_blk_col_start)
     618              : 
     619              :       ! construct the atomic filter matrix
     620          768 :       ALLOCATE (atomic_filter_mat(nrows_atmatrix, ncols_atmatrix))
     621              :       ! calculate atomic filter matrix only if it is non-zero sized
     622          256 :       IF (nrows_atmatrix > 0 .AND. ncols_atmatrix > 0) THEN
     623              :          CALL fb_fltrmat_build_atomic_fltrmat(atomic_H, &
     624              :                                               atomic_S, &
     625              :                                               fermi_level, &
     626              :                                               filter_temp, &
     627              :                                               atomic_filter_mat, &
     628          256 :                                               tolerance)
     629              :       END IF
     630              : 
     631              :       ! ----------------------------------------------------------------------
     632              :       ! Construct filter matrix blocks and add to the correct locations
     633              :       ! in send_buffer
     634              :       ! ----------------------------------------------------------------------
     635              : 
     636              :       ! preconstruct iatom_global to iatom_in_halo map
     637          512 :       ALLOCATE (ind_in_halo(natoms_global))
     638          256 :       ind_in_halo = 0
     639         2304 :       DO iatom_in_halo = 1, natoms_in_halo
     640         2048 :          iatom_global = halo_atoms(iatom_in_halo)
     641         2304 :          ind_in_halo(iatom_global) = iatom_in_halo
     642              :       END DO
     643              : 
     644              :       ! initialise send buffer
     645          256 :       IF (SIZE(send_buf) > 0) send_buf = 0.0_dp
     646              :       ! assign values
     647          768 :       DO ipe = 1, numprocs
     648          512 :          send_sizes(ipe) = 0
     649         2816 :          DO ipair = 1, send_pair_count(ipe)
     650              :             CALL fb_com_atom_pairs_decode(pairs_send(send_pair_disps(ipe) + ipair), &
     651              :                                           pe, iatom_global, jatom_global, &
     652         2048 :                                           send_encode)
     653         2048 :             iatom_in_halo = ind_in_halo(iatom_global)
     654         2048 :             CPASSERT(iatom_in_halo > 0)
     655         2048 :             jatom_in_halo = ind_in_halo(jatom_global)
     656         2048 :             CPASSERT(jatom_in_halo > 0)
     657         2048 :             atomic_kind => particle_set(jatom_global)%atomic_kind
     658              :             CALL get_atomic_kind(atomic_kind=atomic_kind, &
     659         2048 :                                  kind_number=jkind)
     660         2048 :             nrows_blk = row_block_size_data(iatom_global)
     661         2048 :             ncols_blk = ntfns(jkind)
     662              : 
     663              :             ! do it column-wise one trial function at a time
     664        10240 :             DO itrial = 1, ntfns(jkind)
     665         8192 :                ind = send_disps(ipe) + send_sizes(ipe) + (itrial - 1)*nrows_blk
     666              :                CALL dgemv("N", &
     667              :                           nrows_blk, &
     668              :                           ncols_atmatrix, &
     669              :                           1.0_dp, &
     670              :                           atomic_filter_mat( &
     671              :                           atomic_H_blk_row_start(iatom_in_halo): &
     672              :                           atomic_H_blk_row_start(iatom_in_halo + 1) - 1, &
     673              :                           1:ncols_atmatrix &
     674              :                           ), &
     675              :                           nrows_blk, &
     676              :                           atomic_S( &
     677              :                           1:nrows_atmatrix, &
     678              :                           atomic_S_blk_col_start(jatom_in_halo) + &
     679              :                           tfns(itrial, jkind) - 1 &
     680              :                           ), &
     681              :                           1, &
     682              :                           0.0_dp, &
     683              :                           send_buf(ind + 1:ind + nrows_blk), &
     684     23865344 :                           1)
     685              :             END DO ! itrial
     686         6656 :             send_sizes(ipe) = send_sizes(ipe) + nrows_blk*ncols_blk
     687              :          END DO ! ipair
     688              :       END DO ! ipe
     689              : 
     690          256 :       DEALLOCATE (atomic_H)
     691          256 :       DEALLOCATE (atomic_H_blk_row_start)
     692          256 :       DEALLOCATE (atomic_S)
     693          256 :       DEALLOCATE (atomic_S_blk_row_start)
     694          256 :       DEALLOCATE (atomic_filter_mat)
     695          256 :       DEALLOCATE (ind_in_halo)
     696              : 
     697              :       ! ----------------------------------------------------------------------
     698              :       ! Do communication
     699              :       ! ----------------------------------------------------------------------
     700              : 
     701          256 :       CALL timeset("fb_fltrmat_add_blkcol_mpi", handle_mpi)
     702              : 
     703              :       CALL para_env%alltoall(send_buf, send_sizes, send_disps, &
     704          256 :                              recv_buf, recv_sizes, recv_disps)
     705              : 
     706          256 :       CALL timestop(handle_mpi)
     707              : 
     708          256 :       DEALLOCATE (send_buf)
     709          256 :       DEALLOCATE (send_sizes)
     710          256 :       DEALLOCATE (send_disps)
     711          256 :       DEALLOCATE (send_pair_count)
     712          256 :       DEALLOCATE (send_pair_disps)
     713              : 
     714              :       ! ----------------------------------------------------------------------
     715              :       ! Unpack the recv buffer and add the blocks to correct parts of
     716              :       ! the DBCSR filter matrix
     717              :       ! ----------------------------------------------------------------------
     718              : 
     719          768 :       DO ipe = 1, numprocs
     720          512 :          recv_sizes(ipe) = 0
     721         2816 :          DO ipair = 1, recv_pair_count(ipe)
     722              :             CALL fb_com_atom_pairs_decode(pairs_recv(recv_pair_disps(ipe) + ipair), &
     723              :                                           pe, iatom_global, jatom_global, &
     724         2048 :                                           recv_encode)
     725         2048 :             nrows_blk = row_block_size_data(iatom_global)
     726         2048 :             ncols_blk = col_block_size_data(jatom_global)
     727         2048 :             ind = recv_disps(ipe) + recv_sizes(ipe)
     728         2048 :             vector => recv_buf((ind + 1):(ind + nrows_blk*ncols_blk))
     729              :             CALL dbcsr_put_block(filter_mat, &
     730              :                                  iatom_global, jatom_global, &
     731         6144 :                                  block=RESHAPE(vector, [nrows_blk, ncols_blk]))
     732         4608 :             recv_sizes(ipe) = recv_sizes(ipe) + nrows_blk*ncols_blk
     733              :          END DO ! ipair
     734              :       END DO ! ipe
     735              : 
     736              :       ! cleanup rest of the temporary arrays
     737          256 :       DEALLOCATE (recv_buf)
     738          256 :       DEALLOCATE (recv_sizes)
     739          256 :       DEALLOCATE (recv_pair_count)
     740          256 :       DEALLOCATE (recv_pair_disps)
     741              : 
     742          256 :       CALL fb_com_atom_pairs_release(com_pairs_send)
     743          256 :       CALL fb_com_atom_pairs_release(com_pairs_recv)
     744              : 
     745              :       ! cannot finalise the matrix until all blocks has been added
     746              : 
     747          256 :       CALL timestop(handle)
     748              : 
     749         1792 :    END SUBROUTINE fb_fltrmat_add_blkcol
     750              : 
     751              : ! **************************************************************************************************
     752              : !> \brief Computed blocks in one filter matrix column. This version is used by
     753              : !>        fb_fltrmat_build_2, where MPI communication is done collectively
     754              : !> \param H_mat : DBCSR system KS matrix
     755              : !> \param S_mat : DBCSR system overlap matrix
     756              : !> \param H_mat_data  :  local storage of the relevant H_mat matrix blocks
     757              : !> \param S_mat_data  :  local storage of the relevant S_mat matrix blocks
     758              : !> \param atomic_halo :  the halo that contributes to the blk
     759              : !>                       col of the filter matrix
     760              : !> \param trial_fns   :  trial functions data
     761              : !> \param particle_set : set of all particles in the system
     762              : !> \param fermi_level : the fermi level used for defining the filter
     763              : !>                      function, which is a Fermi-Dirac distribution
     764              : !>                      function
     765              : !> \param filter_temp : the filter temperature used for defining the
     766              : !>                      filter function
     767              : !> \param filter_mat_data : local storage for the the computed filter matrix
     768              : !>                          blocks
     769              : !> \param tolerance : anything less than this is regarded as zero
     770              : !> \author Lianheng Tong (LT) lianheng.tong@kcl.ac.uk
     771              : ! **************************************************************************************************
     772           64 :    SUBROUTINE fb_fltrmat_add_blkcol_2(H_mat, &
     773              :                                       S_mat, &
     774              :                                       H_mat_data, &
     775              :                                       S_mat_data, &
     776              :                                       atomic_halo, &
     777              :                                       trial_fns, &
     778              :                                       particle_set, &
     779              :                                       fermi_level, &
     780              :                                       filter_temp, &
     781              :                                       filter_mat_data, &
     782              :                                       tolerance)
     783              :       TYPE(dbcsr_type), POINTER                          :: H_mat, S_mat
     784              :       TYPE(fb_matrix_data_obj), INTENT(IN)               :: H_mat_data, S_mat_data
     785              :       TYPE(fb_atomic_halo_obj), INTENT(IN)               :: atomic_halo
     786              :       TYPE(fb_trial_fns_obj), INTENT(IN)                 :: trial_fns
     787              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     788              :       REAL(KIND=dp), INTENT(IN)                          :: fermi_level, filter_temp
     789              :       TYPE(fb_matrix_data_obj), INTENT(INOUT)            :: filter_mat_data
     790              :       REAL(KIND=dp), INTENT(IN)                          :: tolerance
     791              : 
     792              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'fb_fltrmat_add_blkcol_2'
     793              : 
     794              :       INTEGER :: handle, iatom_global, iatom_in_halo, itrial, jatom_global, jatom_in_halo, jkind, &
     795              :          natoms_global, natoms_in_halo, ncols_atmatrix, ncols_blk, ncols_blk_max, nrows_atmatrix, &
     796              :          nrows_blk, nrows_blk_max
     797           64 :       INTEGER, ALLOCATABLE, DIMENSION(:) :: atomic_H_blk_col_start, atomic_H_blk_row_start, &
     798           64 :          atomic_S_blk_col_start, atomic_S_blk_row_start, col_block_size_data
     799           64 :       INTEGER, DIMENSION(:), POINTER                     :: halo_atoms, ntfns, row_block_size_data
     800           64 :       INTEGER, DIMENSION(:, :), POINTER                  :: tfns
     801              :       LOGICAL                                            :: check_ok
     802           64 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: atomic_filter_mat, atomic_H, atomic_S, &
     803           64 :                                                             mat_blk
     804              :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
     805              : 
     806           64 :       CALL timeset(routineN, handle)
     807              : 
     808           64 :       NULLIFY (atomic_kind, halo_atoms, ntfns, row_block_size_data, tfns)
     809              : 
     810           64 :       check_ok = fb_matrix_data_has_data(H_mat_data)
     811           64 :       CPASSERT(check_ok)
     812           64 :       check_ok = fb_matrix_data_has_data(S_mat_data)
     813           64 :       CPASSERT(check_ok)
     814              : 
     815              :       ! obtain trial function information
     816              :       CALL fb_trial_fns_get(trial_fns=trial_fns, &
     817              :                             nfunctions=ntfns, &
     818           64 :                             functions=tfns)
     819              : 
     820              :       ! obtain row and col block size data for filter matrix
     821           64 :       CALL dbcsr_get_info(H_mat, row_blk_size=row_block_size_data)
     822           64 :       natoms_global = SIZE(particle_set)
     823          192 :       ALLOCATE (col_block_size_data(natoms_global))
     824          576 :       DO jatom_global = 1, natoms_global
     825          512 :          atomic_kind => particle_set(jatom_global)%atomic_kind
     826          512 :          CALL get_atomic_kind(atomic_kind=atomic_kind, kind_number=jkind)
     827          576 :          col_block_size_data(jatom_global) = ntfns(jkind)
     828              :       END DO
     829              : 
     830              :       ! ----------------------------------------------------------------------
     831              :       ! Construct atomic filter matrix for this atomic_halo
     832              :       ! ----------------------------------------------------------------------
     833              : 
     834              :       CALL fb_atomic_halo_get(atomic_halo=atomic_halo, &
     835              :                               natoms=natoms_in_halo, &
     836           64 :                               halo_atoms=halo_atoms)
     837              : 
     838              :       ! construct atomic matrix for H for atomic_halo
     839              :       ALLOCATE (atomic_H_blk_row_start(natoms_in_halo + 1), &
     840          256 :                 atomic_H_blk_col_start(natoms_in_halo + 1))
     841              :       CALL fb_atmatrix_calc_size(H_mat, &
     842              :                                  atomic_halo, &
     843              :                                  nrows_atmatrix, &
     844              :                                  ncols_atmatrix, &
     845              :                                  atomic_H_blk_row_start, &
     846           64 :                                  atomic_H_blk_col_start)
     847          256 :       ALLOCATE (atomic_H(nrows_atmatrix, ncols_atmatrix))
     848              :       CALL fb_atmatrix_construct_2(H_mat_data, &
     849              :                                    atomic_halo, &
     850              :                                    atomic_H, &
     851              :                                    atomic_H_blk_row_start, &
     852           64 :                                    atomic_H_blk_col_start)
     853              : 
     854              :       ! construct atomic matrix for S for atomic_halo
     855              :       ALLOCATE (atomic_S_blk_row_start(natoms_in_halo + 1), &
     856          192 :                 atomic_S_blk_col_start(natoms_in_halo + 1))
     857              :       CALL fb_atmatrix_calc_size(S_mat, &
     858              :                                  atomic_halo, &
     859              :                                  nrows_atmatrix, &
     860              :                                  ncols_atmatrix, &
     861              :                                  atomic_S_blk_row_start, &
     862           64 :                                  atomic_S_blk_col_start)
     863          256 :       ALLOCATE (atomic_S(nrows_atmatrix, ncols_atmatrix))
     864              :       CALL fb_atmatrix_construct_2(S_mat_data, &
     865              :                                    atomic_halo, &
     866              :                                    atomic_S, &
     867              :                                    atomic_S_blk_row_start, &
     868           64 :                                    atomic_S_blk_col_start)
     869              : 
     870              :       ! construct the atomic filter matrix
     871          192 :       ALLOCATE (atomic_filter_mat(nrows_atmatrix, ncols_atmatrix))
     872              :       ! calculate atomic filter matrix only if it is non-zero sized
     873           64 :       IF (nrows_atmatrix > 0 .AND. ncols_atmatrix > 0) THEN
     874              :          CALL fb_fltrmat_build_atomic_fltrmat(atomic_H, &
     875              :                                               atomic_S, &
     876              :                                               fermi_level, &
     877              :                                               filter_temp, &
     878              :                                               atomic_filter_mat, &
     879           64 :                                               tolerance)
     880              :       END IF
     881              : 
     882              :       ! ----------------------------------------------------------------------
     883              :       ! Construct filter matrix block and add to filter_mat_data
     884              :       ! ----------------------------------------------------------------------
     885              : 
     886              :       CALL fb_atomic_halo_get(atomic_halo=atomic_halo, &
     887              :                               owner_atom=jatom_global, &
     888           64 :                               owner_id_in_halo=jatom_in_halo)
     889          576 :       nrows_blk_max = MAXVAL(row_block_size_data)
     890          128 :       ncols_blk_max = MAXVAL(ntfns)
     891          256 :       ALLOCATE (mat_blk(nrows_blk_max, ncols_blk_max))
     892           64 :       mat_blk(:, :) = 0.0_dp
     893          576 :       DO iatom_in_halo = 1, natoms_in_halo
     894          512 :          iatom_global = halo_atoms(iatom_in_halo)
     895          512 :          atomic_kind => particle_set(jatom_global)%atomic_kind
     896              :          CALL get_atomic_kind(atomic_kind=atomic_kind, &
     897          512 :                               kind_number=jkind)
     898          512 :          nrows_blk = row_block_size_data(iatom_global)
     899          512 :          ncols_blk = ntfns(jkind)
     900              : 
     901              :          ! do it column-wise one trial function at a time
     902         2560 :          DO itrial = 1, ntfns(jkind)
     903              :             CALL dgemv("N", &
     904              :                        nrows_blk, &
     905              :                        ncols_atmatrix, &
     906              :                        1.0_dp, &
     907              :                        atomic_filter_mat( &
     908              :                        atomic_H_blk_row_start(iatom_in_halo): &
     909              :                        atomic_H_blk_row_start(iatom_in_halo + 1) - 1, &
     910              :                        1:ncols_atmatrix &
     911              :                        ), &
     912              :                        nrows_blk, &
     913              :                        atomic_S( &
     914              :                        1:nrows_atmatrix, &
     915              :                        atomic_S_blk_col_start(jatom_in_halo) + &
     916              :                        tfns(itrial, jkind) - 1 &
     917              :                        ), &
     918              :                        1, &
     919              :                        0.0_dp, &
     920              :                        mat_blk( &
     921              :                        1:nrows_blk, &
     922              :                        itrial), &
     923      5966336 :                        1)
     924              :          END DO ! itrial
     925              :          CALL fb_matrix_data_add(filter_mat_data, &
     926              :                                  iatom_global, &
     927              :                                  jatom_global, &
     928         1088 :                                  mat_blk(1:nrows_blk, 1:ncols_blk))
     929              :          ! DEALLOCATE(mat_blk, STAT=stat)
     930              :          ! CPPostcondition(stat==0, cp_failure_level, routineP,failure)
     931              :       END DO ! iatom_in_halo
     932           64 :       DEALLOCATE (mat_blk)
     933              : 
     934              :       ! clean up
     935           64 :       DEALLOCATE (atomic_H)
     936           64 :       DEALLOCATE (atomic_H_blk_row_start)
     937           64 :       DEALLOCATE (atomic_S)
     938           64 :       DEALLOCATE (atomic_S_blk_row_start)
     939           64 :       DEALLOCATE (atomic_filter_mat)
     940              : 
     941           64 :       CALL timestop(handle)
     942              : 
     943          384 :    END SUBROUTINE fb_fltrmat_add_blkcol_2
     944              : 
     945              : ! **************************************************************************************************
     946              : !> \brief generate the list of blocks (atom pairs) to be sent and received
     947              : !>        in order to construct the filter matrix for each atomic halo.
     948              : !>        This version is for use with fb_fltrmat_build, where MPI
     949              : !>        communications are done at each step
     950              : !> \param filter_mat : DBCSR formatted filter matrix
     951              : !> \param atomic_halo :  the halo that contributes to a blk
     952              : !>                       col of the filter matrix
     953              : !> \param para_env : cp2k parallel environment
     954              : !> \param atom_pairs_send : list of blocks to be sent
     955              : !> \param atom_pairs_recv : list of blocks to be received
     956              : !> \author Lianheng Tong (LT) lianheng.tong@kcl.ac.uk
     957              : ! **************************************************************************************************
     958          256 :    SUBROUTINE fb_fltrmat_generate_com_pairs(filter_mat, &
     959              :                                             atomic_halo, &
     960              :                                             para_env, &
     961              :                                             atom_pairs_send, &
     962              :                                             atom_pairs_recv)
     963              :       TYPE(dbcsr_type), POINTER                          :: filter_mat
     964              :       TYPE(fb_atomic_halo_obj), INTENT(IN)               :: atomic_halo
     965              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     966              :       TYPE(fb_com_atom_pairs_obj), INTENT(INOUT)         :: atom_pairs_send, atom_pairs_recv
     967              : 
     968              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'fb_fltrmat_generate_com_pairs'
     969              : 
     970              :       INTEGER                                            :: dest, handle, iatom_global, &
     971              :                                                             iatom_in_halo, itask, jatom_global, &
     972              :                                                             natoms_in_halo, nblkrows_total, &
     973              :                                                             ntasks_send
     974          256 :       INTEGER(KIND=int_8), DIMENSION(:, :), POINTER      :: tasks_send
     975          256 :       INTEGER, DIMENSION(:), POINTER                     :: halo_atoms
     976              :       TYPE(fb_com_tasks_obj)                             :: com_tasks_recv, com_tasks_send
     977              : 
     978          256 :       CALL timeset(routineN, handle)
     979              : 
     980          256 :       NULLIFY (tasks_send)
     981          256 :       CALL fb_com_tasks_nullify(com_tasks_send)
     982          256 :       CALL fb_com_tasks_nullify(com_tasks_recv)
     983              : 
     984              :       ! initialise atom_pairs_send and atom_pairs_recv
     985          256 :       IF (fb_com_atom_pairs_has_data(atom_pairs_send)) THEN
     986          256 :          CALL fb_com_atom_pairs_init(atom_pairs_send)
     987              :       ELSE
     988            0 :          CALL fb_com_atom_pairs_create(atom_pairs_send)
     989              :       END IF
     990          256 :       IF (fb_com_atom_pairs_has_data(atom_pairs_recv)) THEN
     991          256 :          CALL fb_com_atom_pairs_init(atom_pairs_recv)
     992              :       ELSE
     993            0 :          CALL fb_com_atom_pairs_create(atom_pairs_recv)
     994              :       END IF
     995              : 
     996              :       ! The total number of filter matrix blocks each processor is going
     997              :       ! to construct equals to the total number of halo atoms in all of
     998              :       ! the atomic halos local to the processor. The number of send
     999              :       ! tasks will not exceed this. We do one halo (col) at a time, and
    1000              :       ! each call of this subroutine will only work on one filter matrix
    1001              :       ! col corresponding to atomic_halo.
    1002              : 
    1003              :       ! The col atom block index for each filter matrix block are the
    1004              :       ! owner atom of each halo. The row atom block index for each
    1005              :       ! filter matrix block corresponding to each col are the halo atoms
    1006              :       ! of the corresponding halos. Filter matrix is non-symmetric: it
    1007              :       ! is non-square, because the blocks themselves are non-square
    1008              : 
    1009              :       CALL fb_atomic_halo_get(atomic_halo=atomic_halo, &
    1010              :                               owner_atom=jatom_global, &
    1011              :                               natoms=natoms_in_halo, &
    1012          256 :                               halo_atoms=halo_atoms)
    1013          256 :       ntasks_send = natoms_in_halo
    1014              : 
    1015              :       ! allocate send tasks
    1016          768 :       ALLOCATE (tasks_send(TASK_N_RECORDS, ntasks_send))
    1017              : 
    1018              :       ! Get the total number of atoms, this can be obtained from the
    1019              :       ! total number of block rows in the DBCSR filter matrix.  We
    1020              :       ! assumes that before calling this subroutine, the filter_mat has
    1021              :       ! already been created and initialised: i.e. using
    1022              :       ! dbcsr_create_new. Even if the matrix is at the moment empty,
    1023              :       ! the attribute nblkrows_total is already assigned from the dbcsr
    1024              :       ! distribution data
    1025              :       CALL dbcsr_get_info(filter_mat, &
    1026          256 :                           nblkrows_total=nblkrows_total)
    1027              : 
    1028              :       ! source is always the local processor
    1029              :       ASSOCIATE (src => para_env%mepos)
    1030              :          ! construct send tasks
    1031          256 :          itask = 1
    1032         2304 :          DO iatom_in_halo = 1, natoms_in_halo
    1033         2048 :             iatom_global = halo_atoms(iatom_in_halo)
    1034              :             ! find where the constructed block of filter matrix belongs to
    1035              :             CALL dbcsr_get_stored_coordinates(filter_mat, &
    1036              :                                               iatom_global, &
    1037              :                                               jatom_global, &
    1038         2048 :                                               processor=dest)
    1039              :             ! create the send tasks
    1040         2048 :             tasks_send(TASK_DEST, itask) = dest
    1041         2048 :             tasks_send(TASK_SRC, itask) = src
    1042              :             CALL fb_com_tasks_encode_pair(tasks_send(TASK_PAIR, itask), &
    1043              :                                           iatom_global, jatom_global, &
    1044         2048 :                                           nblkrows_total)
    1045              :             ! calculation of cost not implemented at the moment
    1046         2048 :             tasks_send(TASK_COST, itask) = 0
    1047         4352 :             itask = itask + 1
    1048              :          END DO ! iatom_in_halo
    1049              :       END ASSOCIATE
    1050              : 
    1051          256 :       CALL fb_com_tasks_create(com_tasks_recv)
    1052          256 :       CALL fb_com_tasks_create(com_tasks_send)
    1053              : 
    1054              :       CALL fb_com_tasks_set(com_tasks=com_tasks_send, &
    1055              :                             task_dim=TASK_N_RECORDS, &
    1056              :                             ntasks=ntasks_send, &
    1057              :                             nencode=nblkrows_total, &
    1058          256 :                             tasks=tasks_send)
    1059              : 
    1060              :       ! generate the recv task list (tasks_recv) from the send task list
    1061              :       CALL fb_com_tasks_transpose_dest_src(com_tasks_recv, "<", com_tasks_send, &
    1062          256 :                                            para_env)
    1063              : 
    1064              :       ! task lists are now complete, now construct the atom_pairs_send
    1065              :       ! and atom_pairs_recv from the tasks lists
    1066              :       CALL fb_com_tasks_build_atom_pairs(com_tasks=com_tasks_send, &
    1067              :                                          atom_pairs=atom_pairs_send, &
    1068              :                                          natoms_encode=nblkrows_total, &
    1069          256 :                                          send_or_recv="send")
    1070              :       CALL fb_com_tasks_build_atom_pairs(com_tasks=com_tasks_recv, &
    1071              :                                          atom_pairs=atom_pairs_recv, &
    1072              :                                          natoms_encode=nblkrows_total, &
    1073          256 :                                          send_or_recv="recv")
    1074              : 
    1075              :       ! cleanup
    1076          256 :       CALL fb_com_tasks_release(com_tasks_recv)
    1077          256 :       CALL fb_com_tasks_release(com_tasks_send)
    1078              : 
    1079          256 :       CALL timestop(handle)
    1080              : 
    1081          768 :    END SUBROUTINE fb_fltrmat_generate_com_pairs
    1082              : 
    1083              : ! **************************************************************************************************
    1084              : !> \brief generate the list of blocks (atom pairs) to be sent and received
    1085              : !>        in order to construct the filter matrix for each atomic halo.
    1086              : !>        This version is for use with fb_fltrmat_build_2, where MPI
    1087              : !>        communications are done collectively.
    1088              : !> \param filter_mat  : DBCSR formatted filter matrix
    1089              : !> \param atomic_halos : set of all local atomic halos contributing to the
    1090              : !>                       filter matrix
    1091              : !> \param para_env : cp2k parallel environment
    1092              : !> \param atom_pairs_send : list of blocks to be sent
    1093              : !> \param atom_pairs_recv : list of blocks to be received
    1094              : !> \author Lianheng Tong (LT) lianheng.tong@kcl.ac.uk
    1095              : ! **************************************************************************************************
    1096           16 :    SUBROUTINE fb_fltrmat_generate_com_pairs_2(filter_mat, &
    1097              :                                               atomic_halos, &
    1098              :                                               para_env, &
    1099              :                                               atom_pairs_send, &
    1100              :                                               atom_pairs_recv)
    1101              :       TYPE(dbcsr_type), POINTER                          :: filter_mat
    1102              :       TYPE(fb_atomic_halo_list_obj), INTENT(IN)          :: atomic_halos
    1103              :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1104              :       TYPE(fb_com_atom_pairs_obj), INTENT(INOUT)         :: atom_pairs_send, atom_pairs_recv
    1105              : 
    1106              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'fb_fltrmat_generate_com_pairs_2'
    1107              : 
    1108              :       INTEGER :: dest, handle, iatom_global, iatom_in_halo, iatom_stored, ihalo, itask, &
    1109              :          jatom_global, jatom_stored, natoms_in_halo, nblkrows_total, nhalos, ntasks_send
    1110           16 :       INTEGER(KIND=int_8), DIMENSION(:, :), POINTER      :: tasks_send
    1111           16 :       INTEGER, DIMENSION(:), POINTER                     :: halo_atoms
    1112              :       LOGICAL                                            :: transpose
    1113           16 :       TYPE(fb_atomic_halo_obj), DIMENSION(:), POINTER    :: halos
    1114              :       TYPE(fb_com_tasks_obj)                             :: com_tasks_recv, com_tasks_send
    1115              : 
    1116           16 :       CALL timeset(routineN, handle)
    1117              : 
    1118           16 :       NULLIFY (tasks_send)
    1119           16 :       CALL fb_com_tasks_nullify(com_tasks_send)
    1120           16 :       CALL fb_com_tasks_nullify(com_tasks_recv)
    1121              : 
    1122              :       ! initialise atom_pairs_send and atom_pairs_recv
    1123           16 :       IF (fb_com_atom_pairs_has_data(atom_pairs_send)) THEN
    1124           16 :          CALL fb_com_atom_pairs_init(atom_pairs_send)
    1125              :       ELSE
    1126            0 :          CALL fb_com_atom_pairs_create(atom_pairs_send)
    1127              :       END IF
    1128           16 :       IF (fb_com_atom_pairs_has_data(atom_pairs_recv)) THEN
    1129           16 :          CALL fb_com_atom_pairs_init(atom_pairs_recv)
    1130              :       ELSE
    1131            0 :          CALL fb_com_atom_pairs_create(atom_pairs_recv)
    1132              :       END IF
    1133              : 
    1134              :       ! The col atom block index for each filter matrix block are the
    1135              :       ! owner atom of each halo. The row atom block index for each
    1136              :       ! filter matrix block corresponding to each col are the halo atoms
    1137              :       ! of the corresponding halos. Filter matrix is non-symmetric: it
    1138              :       ! is non-square, because the blocks themselves are non-square
    1139              : 
    1140              :       CALL fb_atomic_halo_list_get(atomic_halos=atomic_halos, &
    1141              :                                    nhalos=nhalos, &
    1142           16 :                                    halos=halos)
    1143              : 
    1144              :       ! estimate the maximum number of blocks (i.e. atom paris) to send
    1145           16 :       ntasks_send = 0
    1146           80 :       DO ihalo = 1, nhalos
    1147              :          CALL fb_atomic_halo_get(atomic_halo=halos(ihalo), &
    1148           64 :                                  natoms=natoms_in_halo)
    1149           80 :          ntasks_send = ntasks_send + natoms_in_halo
    1150              :       END DO ! ihalo
    1151              : 
    1152              :       ! allocate send tasks
    1153           48 :       ALLOCATE (tasks_send(TASK_N_RECORDS, ntasks_send))
    1154              : 
    1155              :       ! Get the total number of atoms. This can be obtained from the
    1156              :       ! total number of block rows in the DBCSR filter matrix.  We
    1157              :       ! assumes that before calling this subroutine, the filter_mat has
    1158              :       ! already been created and initialised: i.e. using
    1159              :       ! dbcsr_create_new. Even if the matrix is at the moment empty,
    1160              :       ! the attribute nblkrows_total is already assigned from the dbcsr
    1161              :       ! distribution data
    1162              :       CALL dbcsr_get_info(filter_mat, &
    1163           16 :                           nblkrows_total=nblkrows_total)
    1164              : 
    1165              :       ! source is always the local processor
    1166              :       ASSOCIATE (src => para_env%mepos)
    1167              :          ! construct send tasks
    1168           16 :          itask = 1
    1169           80 :          DO ihalo = 1, nhalos
    1170              :             CALL fb_atomic_halo_get(atomic_halo=halos(ihalo), &
    1171              :                                     owner_atom=jatom_global, &
    1172              :                                     natoms=natoms_in_halo, &
    1173           64 :                                     halo_atoms=halo_atoms)
    1174          592 :             DO iatom_in_halo = 1, natoms_in_halo
    1175          512 :                iatom_global = halo_atoms(iatom_in_halo)
    1176          512 :                iatom_stored = iatom_global
    1177          512 :                jatom_stored = jatom_global
    1178          512 :                transpose = .FALSE.
    1179              :                ! find where the constructed block of filter matrix belongs to
    1180              :                CALL dbcsr_get_stored_coordinates(filter_mat, &
    1181              :                                                  iatom_stored, &
    1182              :                                                  jatom_stored, &
    1183          512 :                                                  processor=dest)
    1184              :                ! create the send tasks
    1185          512 :                tasks_send(TASK_DEST, itask) = dest
    1186          512 :                tasks_send(TASK_SRC, itask) = src
    1187              :                CALL fb_com_tasks_encode_pair(tasks_send(TASK_PAIR, itask), &
    1188              :                                              iatom_global, jatom_global, &
    1189          512 :                                              nblkrows_total)
    1190              :                ! calculation of cost not implemented at the moment
    1191          512 :                tasks_send(TASK_COST, itask) = 0
    1192         1088 :                itask = itask + 1
    1193              :             END DO ! iatom_in_halo
    1194              :          END DO ! ihalo
    1195              :       END ASSOCIATE
    1196              : 
    1197              :       ! get the actual number of tasks
    1198           16 :       ntasks_send = itask - 1
    1199              : 
    1200           16 :       CALL fb_com_tasks_create(com_tasks_send)
    1201              :       CALL fb_com_tasks_set(com_tasks=com_tasks_send, &
    1202              :                             task_dim=TASK_N_RECORDS, &
    1203              :                             ntasks=ntasks_send, &
    1204              :                             nencode=nblkrows_total, &
    1205           16 :                             tasks=tasks_send)
    1206              : 
    1207              :       ! generate the recv task list (tasks_recv) from the send task list
    1208           16 :       CALL fb_com_tasks_create(com_tasks_recv)
    1209              :       CALL fb_com_tasks_transpose_dest_src(com_tasks_recv, "<", com_tasks_send, &
    1210           16 :                                            para_env)
    1211              : 
    1212              :       ! task lists are now complete, now construct the atom_pairs_send
    1213              :       ! and atom_pairs_recv from the tasks lists
    1214              :       CALL fb_com_tasks_build_atom_pairs(com_tasks=com_tasks_send, &
    1215              :                                          atom_pairs=atom_pairs_send, &
    1216              :                                          natoms_encode=nblkrows_total, &
    1217           16 :                                          send_or_recv="send")
    1218              :       CALL fb_com_tasks_build_atom_pairs(com_tasks=com_tasks_recv, &
    1219              :                                          atom_pairs=atom_pairs_recv, &
    1220              :                                          natoms_encode=nblkrows_total, &
    1221           16 :                                          send_or_recv="recv")
    1222              : 
    1223              :       ! cleanup
    1224           16 :       CALL fb_com_tasks_release(com_tasks_recv)
    1225           16 :       CALL fb_com_tasks_release(com_tasks_send)
    1226              : 
    1227           16 :       CALL timestop(handle)
    1228              : 
    1229           48 :    END SUBROUTINE fb_fltrmat_generate_com_pairs_2
    1230              : 
    1231              : ! **************************************************************************************************
    1232              : !> \brief Build the atomic filter matrix for each atomic halo
    1233              : !> \param atomic_H : atomic KS matrix
    1234              : !> \param atomic_S : atomic overlap matrix
    1235              : !> \param fermi_level : fermi level used to construct the Fermi-Dirac
    1236              : !>                      filter function
    1237              : !> \param filter_temp : temperature used to construct the Fermi-Dirac
    1238              : !>                      filter function
    1239              : !> \param atomic_filter_mat : the atomic filter matrix
    1240              : !> \param tolerance : anything smaller than tolerance is treated as zero
    1241              : !> \author Lianheng Tong (LT) lianheng.tong@kcl.ac.uk
    1242              : ! **************************************************************************************************
    1243          960 :    SUBROUTINE fb_fltrmat_build_atomic_fltrmat(atomic_H, &
    1244          320 :                                               atomic_S, &
    1245              :                                               fermi_level, &
    1246              :                                               filter_temp, &
    1247          320 :                                               atomic_filter_mat, &
    1248              :                                               tolerance)
    1249              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: atomic_H, atomic_S
    1250              :       REAL(KIND=dp), INTENT(IN)                          :: fermi_level, filter_temp
    1251              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT)        :: atomic_filter_mat
    1252              :       REAL(KIND=dp), INTENT(IN)                          :: tolerance
    1253              : 
    1254              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'fb_fltrmat_build_atomic_fltrmat'
    1255              : 
    1256              :       CHARACTER(LEN=default_string_length)               :: err
    1257              :       INTEGER                                            :: handle, handle_dgemm, handle_dsygv, ii, &
    1258              :                                                             info, jj, mat_dim, work_array_size
    1259              :       LOGICAL                                            :: check_ok
    1260          320 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: eigenvalues, filter_function, work
    1261          320 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: atomic_S_copy, eigenvectors, &
    1262          320 :                                                             filtered_eigenvectors
    1263              : 
    1264          320 :       CALL timeset(routineN, handle)
    1265              : 
    1266              :       ! This subroutine assumes atomic_filter_mat is not zero size, in
    1267              :       ! other words, it really has to be constructed, instead of just
    1268              :       ! being a dummy
    1269              : 
    1270              :       check_ok = SIZE(atomic_filter_mat, 1) > 0 .AND. &
    1271          320 :                  SIZE(atomic_filter_mat, 2) > 0
    1272            0 :       CPASSERT(check_ok)
    1273              : 
    1274              :       ! initialise
    1275      3494720 :       atomic_filter_mat = 0.0_dp
    1276          320 :       mat_dim = SIZE(atomic_H, 1)
    1277              : 
    1278              :       ! diagonalise using LAPACK
    1279          960 :       ALLOCATE (eigenvalues(mat_dim))
    1280              :       ! get optimal work array size
    1281          320 :       ALLOCATE (work(1))
    1282              :       ! dsygv will overwrite part of atomic_H and atomic_S, thus need to copy them
    1283         1280 :       ALLOCATE (atomic_S_copy(SIZE(atomic_S, 1), SIZE(atomic_S, 2)))
    1284      3494720 :       atomic_S_copy(:, :) = atomic_S(:, :)
    1285         1280 :       ALLOCATE (eigenvectors(SIZE(atomic_H, 1), SIZE(atomic_H, 2)))
    1286      3494720 :       eigenvectors(:, :) = atomic_H(:, :)
    1287              : 
    1288          320 :       CALL timeset("fb_atomic_filter_dsygv", handle_dsygv)
    1289              : 
    1290          320 :       info = 0
    1291              :       CALL dsygv(1, 'V', 'U', &
    1292              :                  mat_dim, eigenvectors, mat_dim, &
    1293              :                  atomic_S_copy, mat_dim, eigenvalues, &
    1294          320 :                  work, -1, info)
    1295          320 :       work_array_size = NINT(work(1))
    1296              :       ! now allocate work array
    1297          320 :       DEALLOCATE (work)
    1298          960 :       ALLOCATE (work(work_array_size))
    1299          320 :       work = 0.0_dp
    1300              :       ! do calculation
    1301      3494720 :       atomic_S_copy(:, :) = atomic_S(:, :)
    1302      3494720 :       eigenvectors(:, :) = atomic_H(:, :)
    1303          320 :       info = 0
    1304              :       CALL dsygv(1, 'V', 'U', &
    1305              :                  mat_dim, eigenvectors, mat_dim, &
    1306              :                  atomic_S_copy, mat_dim, eigenvalues, &
    1307          320 :                  work, work_array_size, info)
    1308              :       ! check if diagonalisation is successful
    1309          320 :       IF (info /= 0) THEN
    1310            0 :          WRITE (err, *) "DSYGV failed with non-zero exit status ", info
    1311            0 :          CPABORT(TRIM(err))
    1312              :       END IF
    1313              : 
    1314          320 :       CALL timestop(handle_dsygv)
    1315              : 
    1316          320 :       DEALLOCATE (work)
    1317          320 :       DEALLOCATE (atomic_S_copy)
    1318              : 
    1319              :       ! first get the filter function
    1320          960 :       ALLOCATE (filter_function(mat_dim))
    1321          320 :       filter_function = 0.0_dp
    1322              :       CALL fb_fltrmat_fermi_dirac_mu(filter_function, &
    1323              :                                      eigenvalues, &
    1324              :                                      filter_temp, &
    1325          320 :                                      fermi_level)
    1326          320 :       DEALLOCATE (eigenvalues)
    1327              : 
    1328              :       ! atomic_H has the eigenvectors, construct the version of it
    1329              :       ! filtered through the filter function
    1330         1280 :       ALLOCATE (filtered_eigenvectors(mat_dim, mat_dim))
    1331        33600 :       DO jj = 1, mat_dim
    1332      3494720 :          DO ii = 1, mat_dim
    1333              :             filtered_eigenvectors(ii, jj) = &
    1334      3494400 :                filter_function(jj)*eigenvectors(ii, jj)
    1335              :          END DO ! ii
    1336              :       END DO ! jj
    1337              : 
    1338          320 :       DEALLOCATE (filter_function)
    1339              : 
    1340          320 :       CALL timeset("fb_atomic_filter_dgemm", handle_dgemm)
    1341              : 
    1342              :       ! construct atomic filter matrix
    1343              :       CALL dgemm("N", &
    1344              :                  "T", &
    1345              :                  mat_dim, &
    1346              :                  mat_dim, &
    1347              :                  mat_dim, &
    1348              :                  1.0_dp, &
    1349              :                  filtered_eigenvectors, &
    1350              :                  mat_dim, &
    1351              :                  eigenvectors, &
    1352              :                  mat_dim, &
    1353              :                  0.0_dp, &
    1354              :                  atomic_filter_mat, &
    1355          320 :                  mat_dim)
    1356              : 
    1357          320 :       CALL timestop(handle_dgemm)
    1358              : 
    1359              :       ! remove small negative terms due to numerical error, the filter
    1360              :       ! matrix must not be negative definite
    1361        33600 :       DO jj = 1, SIZE(atomic_filter_mat, 2)
    1362      3494720 :          DO ii = 1, SIZE(atomic_filter_mat, 1)
    1363      3494400 :             IF (ABS(atomic_filter_mat(ii, jj)) < tolerance) THEN
    1364        52688 :                atomic_filter_mat(ii, jj) = 0.0_dp
    1365              :             END IF
    1366              :          END DO
    1367              :       END DO
    1368              : 
    1369          320 :       DEALLOCATE (filtered_eigenvectors)
    1370          320 :       DEALLOCATE (eigenvectors)
    1371              : 
    1372          320 :       CALL timestop(handle)
    1373              : 
    1374          960 :    END SUBROUTINE fb_fltrmat_build_atomic_fltrmat
    1375              : 
    1376              : ! **************************************************************************************************
    1377              : !> \brief get values of Fermi-Dirac distribution based on a given fermi
    1378              : !>        level at a given set of energy eigenvalues
    1379              : !> \param f : the Fermi-Dirac distribution function values
    1380              : !> \param eigenvals : set of energy eigenvalues
    1381              : !> \param T : temperature
    1382              : !> \param mu : the fermi level
    1383              : !> \author Lianheng Tong (LT) lianheng.tong@kcl.ac.uk
    1384              : ! **************************************************************************************************
    1385          320 :    SUBROUTINE fb_fltrmat_fermi_dirac_mu(f, eigenvals, T, mu)
    1386              :       REAL(KIND=dp), DIMENSION(:), INTENT(OUT)           :: f
    1387              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: eigenvals
    1388              :       REAL(KIND=dp), INTENT(IN)                          :: T, mu
    1389              : 
    1390              :       REAL(KIND=dp)                                      :: kTS, ne
    1391              : 
    1392              : ! we want fermi function max at 1, so maxocc = 1 here
    1393              : 
    1394          320 :       CALL SmearOcc(f, ne, kTS, eigenvals, mu, T, 1.0_dp, smear_fermi_dirac)
    1395          320 :    END SUBROUTINE fb_fltrmat_fermi_dirac_mu
    1396              : 
    1397              : ! **************************************************************************************************
    1398              : !> \brief get values of Fermi-Dirac distribution based on a given electron
    1399              : !>        number at a given set of energy eigenvales
    1400              : !> \param f : the Fermi-Dirac distribution function values
    1401              : !> \param eigenvals : set of energy eigenvalues
    1402              : !> \param T : temperature
    1403              : !> \param ne : number of electrons
    1404              : !> \param maxocc : maximum occupancy per orbital
    1405              : !> \author Lianheng Tong (LT) lianheng.tong@kcl.ac.uk
    1406              : ! **************************************************************************************************
    1407            0 :    SUBROUTINE fb_fltrmat_fermi_dirac_ne(f, eigenvals, T, ne, maxocc)
    1408              :       REAL(KIND=dp), DIMENSION(:), INTENT(OUT)           :: f
    1409              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: eigenvals
    1410              :       REAL(KIND=dp), INTENT(IN)                          :: T, ne, maxocc
    1411              : 
    1412              :       REAL(KIND=dp)                                      :: kTS, mu
    1413              : 
    1414              : ! mu is the calculated fermi level
    1415              : ! kTS is the calculated entropic contribution to the energy i.e. -TS
    1416              : ! kTS = kT*[f ln f + (1-f) ln (1-f)]
    1417              : 
    1418            0 :       CALL SmearFixed(f, mu, kTS, eigenvals, ne, T, maxocc, smear_fermi_dirac)
    1419            0 :    END SUBROUTINE fb_fltrmat_fermi_dirac_ne
    1420              : 
    1421              : END MODULE qs_fb_filter_matrix_methods
        

Generated by: LCOV version 2.0-1