LCOV - code coverage report
Current view: top level - src - qs_scf_csr_write.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 96.8 % 216 209
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 5 5

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Functions to print the KS and S matrix in the CSR format to file
      10              : !> \par History
      11              : !>      Started as a copy from the relevant part of qs_scf_post_gpw
      12              : !> \author Fabian Ducry (05.2020)
      13              : ! **************************************************************************************************
      14              : MODULE qs_scf_csr_write
      15              :    USE cp_dbcsr_api,                    ONLY: &
      16              :         dbcsr_convert_dbcsr_to_csr, dbcsr_copy, dbcsr_create, &
      17              :         dbcsr_csr_create_and_convert_complex, dbcsr_csr_create_from_dbcsr, &
      18              :         dbcsr_csr_dbcsr_blkrow_dist, dbcsr_csr_destroy, dbcsr_csr_type, dbcsr_csr_write, &
      19              :         dbcsr_desymmetrize, dbcsr_finalize, dbcsr_get_block_p, dbcsr_has_symmetry, dbcsr_p_type, &
      20              :         dbcsr_put_block, dbcsr_release, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, &
      21              :         dbcsr_type_no_symmetry, dbcsr_type_symmetric
      22              :    USE cp_dbcsr_cp2k_link,              ONLY: cp_dbcsr_alloc_block_from_nbl
      23              :    USE cp_dbcsr_operations,             ONLY: dbcsr_allocate_matrix_set,&
      24              :                                               dbcsr_deallocate_matrix_set
      25              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      26              :                                               cp_logger_get_default_io_unit,&
      27              :                                               cp_logger_type
      28              :    USE cp_output_handling,              ONLY: cp_p_file,&
      29              :                                               cp_print_key_finished_output,&
      30              :                                               cp_print_key_should_output,&
      31              :                                               cp_print_key_unit_nr
      32              :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      33              :                                               section_vals_type,&
      34              :                                               section_vals_val_get
      35              :    USE kinds,                           ONLY: default_path_length,&
      36              :                                               dp
      37              :    USE kpoint_methods,                  ONLY: rskp_transform
      38              :    USE kpoint_types,                    ONLY: get_kpoint_info,&
      39              :                                               kpoint_type
      40              :    USE qs_environment_types,            ONLY: get_qs_env,&
      41              :                                               qs_environment_type
      42              :    USE qs_neighbor_list_types,          ONLY: get_iterator_info,&
      43              :                                               get_neighbor_list_set_p,&
      44              :                                               neighbor_list_iterate,&
      45              :                                               neighbor_list_iterator_create,&
      46              :                                               neighbor_list_iterator_p_type,&
      47              :                                               neighbor_list_iterator_release,&
      48              :                                               neighbor_list_set_p_type
      49              : #include "./base/base_uses.f90"
      50              : 
      51              :    IMPLICIT NONE
      52              :    PRIVATE
      53              : 
      54              :    ! Global parameters
      55              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf_csr_write'
      56              :    PUBLIC :: write_ks_matrix_csr, &
      57              :              write_s_matrix_csr
      58              : 
      59              : ! **************************************************************************************************
      60              : 
      61              : CONTAINS
      62              : 
      63              : !**************************************************************************************************
      64              : !> \brief writing the KS matrix in csr format into a file
      65              : !> \param qs_env qs environment
      66              : !> \param input the input
      67              : !> \par History
      68              : !>       Moved to module qs_scf_csr_write (05.2020)
      69              : !> \author Mohammad Hossein Bani-Hashemian
      70              : ! **************************************************************************************************
      71        18069 :    SUBROUTINE write_ks_matrix_csr(qs_env, input)
      72              :       TYPE(qs_environment_type), POINTER                 :: qs_env
      73              :       TYPE(section_vals_type), POINTER                   :: input
      74              : 
      75              :       CHARACTER(len=*), PARAMETER :: routineN = 'write_ks_matrix_csr'
      76              : 
      77              :       INTEGER                                            :: handle, output_unit
      78              :       LOGICAL                                            :: do_kpoints, do_ks_csr_write, real_space
      79              :       TYPE(cp_logger_type), POINTER                      :: logger
      80        18069 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_ks
      81              :       TYPE(kpoint_type), POINTER                         :: kpoints
      82              :       TYPE(section_vals_type), POINTER                   :: dft_section
      83              : 
      84        18069 :       CALL timeset(routineN, handle)
      85              : 
      86              :       NULLIFY (dft_section)
      87              : 
      88        18069 :       logger => cp_get_default_logger()
      89        18069 :       output_unit = cp_logger_get_default_io_unit(logger)
      90              : 
      91        18069 :       dft_section => section_vals_get_subs_vals(input, "DFT")
      92              :       do_ks_csr_write = BTEST(cp_print_key_should_output(logger%iter_info, dft_section, &
      93        18069 :                                                          "PRINT%KS_CSR_WRITE"), cp_p_file)
      94              : 
      95        18069 :       IF (do_ks_csr_write) THEN
      96           54 :          CALL get_qs_env(qs_env=qs_env, kpoints=kpoints, matrix_ks_kp=matrix_ks, do_kpoints=do_kpoints)
      97              :          CALL section_vals_val_get(dft_section, "PRINT%KS_CSR_WRITE%REAL_SPACE", &
      98           54 :                                    l_val=real_space)
      99              : 
     100           54 :          IF (do_kpoints .AND. .NOT. real_space) THEN
     101              :             CALL write_matrix_kp_csr(mat=matrix_ks, dft_section=dft_section, &
     102           10 :                                      kpoints=kpoints, prefix="KS")
     103              :          ELSE
     104              :             CALL write_matrix_csr(dft_section, mat=matrix_ks, kpoints=kpoints, do_kpoints=do_kpoints, &
     105           44 :                                   prefix="KS")
     106              :          END IF
     107              :       END IF
     108              : 
     109        18069 :       CALL timestop(handle)
     110              : 
     111        18069 :    END SUBROUTINE write_ks_matrix_csr
     112              : 
     113              : !**************************************************************************************************
     114              : !> \brief writing the overlap matrix in csr format into a file
     115              : !> \param qs_env qs environment
     116              : !> \param input the input
     117              : !> \par History
     118              : !>      Moved to module qs_scf_csr_write
     119              : !> \author Mohammad Hossein Bani-Hashemian
     120              : ! **************************************************************************************************
     121        18069 :    SUBROUTINE write_s_matrix_csr(qs_env, input)
     122              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     123              :       TYPE(section_vals_type), POINTER                   :: input
     124              : 
     125              :       CHARACTER(len=*), PARAMETER :: routineN = 'write_s_matrix_csr'
     126              : 
     127              :       INTEGER                                            :: handle, output_unit
     128              :       LOGICAL                                            :: do_kpoints, do_s_csr_write, real_space
     129              :       TYPE(cp_logger_type), POINTER                      :: logger
     130        18069 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_s
     131              :       TYPE(kpoint_type), POINTER                         :: kpoints
     132              :       TYPE(section_vals_type), POINTER                   :: dft_section
     133              : 
     134        18069 :       CALL timeset(routineN, handle)
     135              : 
     136              :       NULLIFY (dft_section)
     137              : 
     138        18069 :       logger => cp_get_default_logger()
     139        18069 :       output_unit = cp_logger_get_default_io_unit(logger)
     140              : 
     141        18069 :       dft_section => section_vals_get_subs_vals(input, "DFT")
     142              :       do_s_csr_write = BTEST(cp_print_key_should_output(logger%iter_info, dft_section, &
     143        18069 :                                                         "PRINT%S_CSR_WRITE"), cp_p_file)
     144              : 
     145        18069 :       IF (do_s_csr_write) THEN
     146           52 :          CALL get_qs_env(qs_env=qs_env, kpoints=kpoints, matrix_s_kp=matrix_s, do_kpoints=do_kpoints)
     147              :          CALL section_vals_val_get(dft_section, "PRINT%S_CSR_WRITE%REAL_SPACE", &
     148           52 :                                    l_val=real_space)
     149              : 
     150           52 :          IF (do_kpoints .AND. .NOT. real_space) THEN
     151              :             CALL write_matrix_kp_csr(mat=matrix_s, dft_section=dft_section, &
     152           10 :                                      kpoints=kpoints, prefix="S")
     153              :          ELSE
     154              :             CALL write_matrix_csr(dft_section, mat=matrix_s, kpoints=kpoints, do_kpoints=do_kpoints, &
     155           42 :                                   prefix="S")
     156              :          END IF
     157              :       END IF
     158              : 
     159        18069 :       CALL timestop(handle)
     160              : 
     161        18069 :    END SUBROUTINE write_s_matrix_csr
     162              : 
     163              : ! **************************************************************************************************
     164              : !> \brief helper function to print the real space representation of KS or S matrix to file
     165              : !> \param dft_section the dft_section
     166              : !> \param mat Hamiltonian or overlap matrix
     167              : !> \param kpoints Kpoint environment
     168              : !> \param prefix string to distinguish between KS and S matrix
     169              : !> \param do_kpoints Whether it is a gamma-point or k-point calculation
     170              : !> \par History
     171              : !>       Moved most of the code from write_ks_matrix_csr and write_s_matrix_csr
     172              : !>       Removed the code for printing k-point dependent matrices and added
     173              : !>              printing of the real space representation
     174              : ! **************************************************************************************************
     175           86 :    SUBROUTINE write_matrix_csr(dft_section, mat, kpoints, prefix, do_kpoints)
     176              :       TYPE(section_vals_type), INTENT(IN), POINTER       :: dft_section
     177              :       TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(IN), &
     178              :          POINTER                                         :: mat
     179              :       TYPE(kpoint_type), INTENT(IN), POINTER             :: kpoints
     180              :       CHARACTER(*), INTENT(in)                           :: prefix
     181              :       LOGICAL, INTENT(IN)                                :: do_kpoints
     182              : 
     183              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'write_matrix_csr'
     184              : 
     185              :       CHARACTER(LEN=default_path_length)                 :: file_name, fileformat, subs_string
     186              :       INTEGER                                            :: handle, ic, ispin, ncell, nspin, &
     187              :                                                             output_unit, unit_nr
     188           86 :       INTEGER, ALLOCATABLE, DIMENSION(:, :)              :: index_to_cell
     189           86 :       INTEGER, ALLOCATABLE, DIMENSION(:, :, :)           :: cell_to_index
     190           86 :       INTEGER, DIMENSION(:, :), POINTER                  :: i2c
     191           86 :       INTEGER, DIMENSION(:, :, :), POINTER               :: c2i
     192              :       LOGICAL                                            :: bin, do_symmetric, uptr
     193              :       REAL(KIND=dp)                                      :: thld
     194              :       TYPE(cp_logger_type), POINTER                      :: logger
     195              :       TYPE(dbcsr_csr_type)                               :: mat_csr
     196           86 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: mat_nosym
     197              :       TYPE(dbcsr_type)                                   :: matrix_nosym
     198              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     199           86 :          POINTER                                         :: sab_nl
     200              : 
     201           86 :       CALL timeset(routineN, handle)
     202              : 
     203           86 :       logger => cp_get_default_logger()
     204           86 :       output_unit = cp_logger_get_default_io_unit(logger)
     205              : 
     206           86 :       subs_string = "PRINT%"//prefix//"_CSR_WRITE"
     207              : 
     208           86 :       CALL section_vals_val_get(dft_section, subs_string//"%THRESHOLD", r_val=thld)
     209           86 :       CALL section_vals_val_get(dft_section, subs_string//"%UPPER_TRIANGULAR", l_val=uptr)
     210           86 :       CALL section_vals_val_get(dft_section, subs_string//"%BINARY", l_val=bin)
     211              : 
     212           86 :       IF (bin) THEN
     213            2 :          fileformat = "UNFORMATTED"
     214              :       ELSE
     215           84 :          fileformat = "FORMATTED"
     216              :       END IF
     217              : 
     218           86 :       nspin = SIZE(mat, 1)
     219           86 :       ncell = SIZE(mat, 2)
     220              : 
     221           86 :       IF (do_kpoints) THEN
     222              : 
     223            2 :          i2c => kpoints%index_to_cell
     224            2 :          c2i => kpoints%cell_to_index
     225              : 
     226            2 :          NULLIFY (sab_nl)
     227            2 :          CALL get_kpoint_info(kpoints, sab_nl=sab_nl)
     228            2 :          CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
     229              : 
     230              :          ! desymmetrize the KS or S matrices if necessary
     231            2 :          IF (do_symmetric) THEN
     232              :             CALL desymmetrize_rs_matrix(mat, mat_nosym, cell_to_index, index_to_cell, kpoints)
     233            2 :             ncell = SIZE(index_to_cell, 2) ! update the number of cells
     234              :          ELSE
     235              :             ALLOCATE (cell_to_index(LBOUND(c2i, 1):UBOUND(c2i, 1), &
     236              :                                     LBOUND(c2i, 2):UBOUND(c2i, 2), &
     237            0 :                                     LBOUND(c2i, 3):UBOUND(c2i, 3)))
     238              :             cell_to_index(LBOUND(c2i, 1):UBOUND(c2i, 1), &
     239              :                           LBOUND(c2i, 2):UBOUND(c2i, 2), &
     240            0 :                           LBOUND(c2i, 3):UBOUND(c2i, 3)) = c2i
     241              : 
     242            0 :             ALLOCATE (index_to_cell(3, ncell))
     243            0 :             index_to_cell(1:3, 1:ncell) = i2c
     244              : 
     245            0 :             mat_nosym => mat
     246              :          END IF
     247              : 
     248              :          ! print the index to cell mapping to the output
     249            2 :          IF (output_unit > 0) THEN
     250            1 :             WRITE (output_unit, "(/,A15,T15,I4,A)") prefix//" CSR write| ", &
     251            2 :                ncell, " periodic images"
     252            1 :             WRITE (output_unit, "(T7,A,T17,A,T24,A,T31,A)") "Number", "X", "Y", "Z"
     253           28 :             DO ic = 1, ncell
     254           28 :                WRITE (output_unit, "(T8,I3,T15,I3,T22,I3,T29,I3)") ic, index_to_cell(:, ic)
     255              :             END DO
     256              :          END IF
     257              :       END IF
     258              : 
     259              :       ! write the csr file(s)
     260          172 :       DO ispin = 1, nspin
     261          310 :          DO ic = 1, ncell
     262          138 :             IF (do_kpoints) THEN
     263           54 :                CALL dbcsr_copy(matrix_nosym, mat_nosym(ispin, ic)%matrix)
     264           54 :                WRITE (file_name, '(2(A,I0))') prefix//"_SPIN_", ispin, "_R_", ic
     265              :             ELSE
     266           84 :                IF (dbcsr_has_symmetry(mat(ispin, ic)%matrix)) THEN
     267           84 :                   CALL dbcsr_desymmetrize(mat(ispin, ic)%matrix, matrix_nosym)
     268              :                ELSE
     269            0 :                   CALL dbcsr_copy(matrix_nosym, mat(ispin, ic)%matrix)
     270              :                END IF
     271           84 :                WRITE (file_name, '(A,I0)') prefix//"_SPIN_", ispin
     272              :             END IF
     273              :             ! Convert dbcsr to csr
     274              :             CALL dbcsr_csr_create_from_dbcsr(matrix_nosym, &
     275          138 :                                              mat_csr, dbcsr_csr_dbcsr_blkrow_dist)
     276          138 :             CALL dbcsr_convert_dbcsr_to_csr(matrix_nosym, mat_csr)
     277              :             ! Write to file
     278              :             unit_nr = cp_print_key_unit_nr(logger, dft_section, subs_string, &
     279              :                                            extension=".csr", middle_name=TRIM(file_name), &
     280          138 :                                            file_status="REPLACE", file_form=fileformat)
     281          138 :             CALL dbcsr_csr_write(mat_csr, unit_nr, upper_triangle=uptr, threshold=thld, binary=bin)
     282              : 
     283          138 :             CALL cp_print_key_finished_output(unit_nr, logger, dft_section, subs_string)
     284          138 :             CALL dbcsr_csr_destroy(mat_csr)
     285          224 :             CALL dbcsr_release(matrix_nosym)
     286              :          END DO
     287              :       END DO
     288              : 
     289              :       ! clean up
     290           86 :       IF (do_kpoints) THEN
     291            2 :          DEALLOCATE (cell_to_index, index_to_cell)
     292            2 :          IF (do_symmetric) THEN
     293            4 :             DO ispin = 1, nspin
     294           58 :                DO ic = 1, ncell
     295           56 :                   CALL dbcsr_release(mat_nosym(ispin, ic)%matrix)
     296              :                END DO
     297              :             END DO
     298            2 :             CALL dbcsr_deallocate_matrix_set(mat_nosym)
     299              :          END IF
     300              :       END IF
     301           86 :       CALL timestop(handle)
     302              : 
     303          172 :    END SUBROUTINE write_matrix_csr
     304              : 
     305              : ! **************************************************************************************************
     306              : !> \brief helper function to print the k-dependent KS or S matrix to file
     307              : !> \param mat Hamiltonian or overlap matrix for k-point calculations
     308              : !> \param dft_section the dft_section
     309              : !> \param kpoints Kpoint environment
     310              : !> \param prefix string to distinguish between KS and S matrix
     311              : !> \par History
     312              : !>       Moved the code from write_matrix_csr to write_matrix_kp_csr
     313              : !> \author Fabian Ducry
     314              : ! **************************************************************************************************
     315           20 :    SUBROUTINE write_matrix_kp_csr(mat, dft_section, kpoints, prefix)
     316              :       TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(IN), &
     317              :          POINTER                                         :: mat
     318              :       TYPE(section_vals_type), INTENT(IN), POINTER       :: dft_section
     319              :       TYPE(kpoint_type), INTENT(IN), POINTER             :: kpoints
     320              :       CHARACTER(*), INTENT(in)                           :: prefix
     321              : 
     322              :       CHARACTER(len=*), PARAMETER :: routineN = 'write_matrix_kp_csr'
     323              : 
     324              :       CHARACTER(LEN=default_path_length)                 :: file_name, fileformat, subs_string
     325              :       INTEGER                                            :: handle, igroup, ik, ikp, ispin, kplocal, &
     326              :                                                             nkp_groups, output_unit, unit_nr
     327              :       INTEGER, DIMENSION(2)                              :: kp_range
     328           20 :       INTEGER, DIMENSION(:, :), POINTER                  :: kp_dist
     329              :       LOGICAL                                            :: bin, uptr, use_real_wfn
     330              :       REAL(KIND=dp)                                      :: thld
     331           20 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: xkp
     332              :       TYPE(cp_logger_type), POINTER                      :: logger
     333              :       TYPE(dbcsr_csr_type)                               :: mat_csr
     334              :       TYPE(dbcsr_type)                                   :: matrix_nosym
     335              :       TYPE(dbcsr_type), POINTER                          :: imatrix, imatrix_nosym, rmatrix, &
     336              :                                                             rmatrix_nosym
     337              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     338           20 :          POINTER                                         :: sab_nl
     339              : 
     340           20 :       CALL timeset(routineN, handle)
     341              : 
     342           20 :       logger => cp_get_default_logger()
     343           20 :       output_unit = cp_logger_get_default_io_unit(logger)
     344              : 
     345           20 :       subs_string = "PRINT%"//prefix//"_CSR_WRITE"
     346              : 
     347           20 :       CALL section_vals_val_get(dft_section, subs_string//"%THRESHOLD", r_val=thld)
     348           20 :       CALL section_vals_val_get(dft_section, subs_string//"%UPPER_TRIANGULAR", l_val=uptr)
     349           20 :       CALL section_vals_val_get(dft_section, subs_string//"%BINARY", l_val=bin)
     350              : 
     351           20 :       IF (bin) THEN
     352           20 :          fileformat = "UNFORMATTED"
     353              :       ELSE
     354            0 :          fileformat = "FORMATTED"
     355              :       END IF
     356              : 
     357           20 :       NULLIFY (sab_nl)
     358              : 
     359              :       !  Calculate the Hamiltonian at the k-points
     360              :       CALL get_kpoint_info(kpoints, xkp=xkp, use_real_wfn=use_real_wfn, kp_range=kp_range, &
     361           20 :                            nkp_groups=nkp_groups, kp_dist=kp_dist, sab_nl=sab_nl)
     362              : 
     363           20 :       ALLOCATE (rmatrix)
     364              :       CALL dbcsr_create(rmatrix, template=mat(1, 1)%matrix, &
     365           20 :                         matrix_type=dbcsr_type_symmetric)
     366           20 :       CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_nl)
     367              : 
     368           20 :       IF (.NOT. use_real_wfn) THEN
     369              :          ! Allocate temporary variables
     370           16 :          ALLOCATE (rmatrix_nosym, imatrix, imatrix_nosym)
     371              :          CALL dbcsr_create(rmatrix_nosym, template=mat(1, 1)%matrix, &
     372           16 :                            matrix_type=dbcsr_type_no_symmetry)
     373              :          CALL dbcsr_create(imatrix, template=mat(1, 1)%matrix, &
     374           16 :                            matrix_type=dbcsr_type_antisymmetric)
     375              :          CALL dbcsr_create(imatrix_nosym, template=mat(1, 1)%matrix, &
     376           16 :                            matrix_type=dbcsr_type_no_symmetry)
     377           16 :          CALL cp_dbcsr_alloc_block_from_nbl(rmatrix_nosym, sab_nl)
     378           16 :          CALL cp_dbcsr_alloc_block_from_nbl(imatrix, sab_nl)
     379           16 :          CALL cp_dbcsr_alloc_block_from_nbl(imatrix_nosym, sab_nl)
     380              :       END IF
     381              : 
     382           20 :       kplocal = kp_range(2) - kp_range(1) + 1
     383           56 :       DO ikp = 1, kplocal
     384           92 :          DO ispin = 1, SIZE(mat, 1)
     385          140 :             DO igroup = 1, nkp_groups
     386              :                ! number of current kpoint
     387           68 :                ik = kp_dist(1, igroup) + ikp - 1
     388           68 :                CALL dbcsr_set(rmatrix, 0.0_dp)
     389           68 :                IF (use_real_wfn) THEN
     390              :                   ! FT of the matrix
     391              :                   CALL rskp_transform(rmatrix=rmatrix, rsmat=mat, ispin=ispin, &
     392            4 :                                       xkp=xkp(1:3, ik), cell_to_index=kpoints%cell_to_index, sab_nl=sab_nl)
     393              :                   ! Convert to desymmetrized csr matrix
     394            4 :                   CALL dbcsr_desymmetrize(rmatrix, matrix_nosym)
     395            4 :                   CALL dbcsr_csr_create_from_dbcsr(matrix_nosym, mat_csr, dbcsr_csr_dbcsr_blkrow_dist)
     396            4 :                   CALL dbcsr_convert_dbcsr_to_csr(matrix_nosym, mat_csr)
     397            4 :                   CALL dbcsr_release(matrix_nosym)
     398              :                ELSE
     399              :                   ! FT of the matrix
     400           64 :                   CALL dbcsr_set(imatrix, 0.0_dp)
     401              :                   CALL rskp_transform(rmatrix=rmatrix, cmatrix=imatrix, rsmat=mat, ispin=ispin, &
     402           64 :                                       xkp=xkp(1:3, ik), cell_to_index=kpoints%cell_to_index, sab_nl=sab_nl)
     403              : 
     404              :                   ! Desymmetrize and sum the real and imaginary part into
     405              :                   ! cmatrix
     406           64 :                   CALL dbcsr_desymmetrize(rmatrix, rmatrix_nosym)
     407           64 :                   CALL dbcsr_desymmetrize(imatrix, imatrix_nosym)
     408              :                   ! Convert to csr
     409              :                   CALL dbcsr_csr_create_and_convert_complex(rmatrix=rmatrix_nosym, &
     410              :                                                             imatrix=imatrix_nosym, &
     411              :                                                             csr_mat=mat_csr, &
     412           64 :                                                             dist_format=dbcsr_csr_dbcsr_blkrow_dist)
     413              :                END IF
     414              :                ! Write to file
     415           68 :                WRITE (file_name, '(2(A,I0))') prefix//"_SPIN_", ispin, "_K_", ik
     416              :                unit_nr = cp_print_key_unit_nr(logger, dft_section, subs_string, &
     417              :                                               extension=".csr", middle_name=TRIM(file_name), &
     418           68 :                                               file_status="REPLACE", file_form=fileformat)
     419           68 :                CALL dbcsr_csr_write(mat_csr, unit_nr, upper_triangle=uptr, threshold=thld, binary=bin)
     420              : 
     421           68 :                CALL cp_print_key_finished_output(unit_nr, logger, dft_section, subs_string)
     422              : 
     423          104 :                CALL dbcsr_csr_destroy(mat_csr)
     424              :             END DO
     425              :          END DO
     426              :       END DO
     427           20 :       CALL dbcsr_release(rmatrix)
     428           20 :       DEALLOCATE (rmatrix)
     429           20 :       IF (.NOT. use_real_wfn) THEN
     430           16 :          CALL dbcsr_release(rmatrix_nosym)
     431           16 :          CALL dbcsr_release(imatrix)
     432           16 :          CALL dbcsr_release(imatrix_nosym)
     433           16 :          DEALLOCATE (rmatrix_nosym, imatrix, imatrix_nosym)
     434              :       END IF
     435           20 :       CALL timestop(handle)
     436              : 
     437           20 :    END SUBROUTINE write_matrix_kp_csr
     438              : 
     439              : ! **************************************************************************************************
     440              : !> \brief Desymmetrizes the KS or S matrices which are stored in symmetric !matrices
     441              : !> \param mat Hamiltonian or overlap matrices
     442              : !> \param mat_nosym Desymmetrized Hamiltonian or overlap matrices
     443              : !> \param cell_to_index Mapping of cell indices to linear RS indices
     444              : !> \param index_to_cell Mapping of linear RS indices to cell indices
     445              : !> \param kpoints Kpoint environment
     446              : !> \author Fabian Ducry
     447              : ! **************************************************************************************************
     448            2 :    SUBROUTINE desymmetrize_rs_matrix(mat, mat_nosym, cell_to_index, index_to_cell, kpoints)
     449              :       TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(IN), &
     450              :          POINTER                                         :: mat
     451              :       TYPE(dbcsr_p_type), DIMENSION(:, :), &
     452              :          INTENT(INOUT), POINTER                          :: mat_nosym
     453              :       INTEGER, ALLOCATABLE, DIMENSION(:, :, :), &
     454              :          INTENT(OUT)                                     :: cell_to_index
     455              :       INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(OUT) :: index_to_cell
     456              :       TYPE(kpoint_type), INTENT(IN), POINTER             :: kpoints
     457              : 
     458              :       CHARACTER(len=*), PARAMETER :: routineN = 'desymmetrize_rs_matrix'
     459              : 
     460              :       INTEGER                                            :: handle, iatom, ic, icn, icol, irow, &
     461              :                                                             ispin, jatom, ncell, nomirror, nspin, &
     462              :                                                             nx, ny, nz
     463              :       INTEGER, DIMENSION(3)                              :: cell
     464            2 :       INTEGER, DIMENSION(:, :), POINTER                  :: i2c
     465            2 :       INTEGER, DIMENSION(:, :, :), POINTER               :: c2i
     466              :       LOGICAL                                            :: found, lwtr
     467            2 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: block
     468              :       TYPE(neighbor_list_iterator_p_type), &
     469            2 :          DIMENSION(:), POINTER                           :: nl_iterator
     470              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     471            2 :          POINTER                                         :: sab_nl
     472              : 
     473            2 :       CALL timeset(routineN, handle)
     474              : 
     475            2 :       i2c => kpoints%index_to_cell
     476            2 :       c2i => kpoints%cell_to_index
     477              : 
     478            2 :       ncell = SIZE(i2c, 2)
     479            2 :       nspin = SIZE(mat, 1)
     480              : 
     481            2 :       nx = MAX(ABS(LBOUND(c2i, 1)), ABS(UBOUND(c2i, 1)))
     482            2 :       ny = MAX(ABS(LBOUND(c2i, 2)), ABS(UBOUND(c2i, 3)))
     483            2 :       nz = MAX(ABS(LBOUND(c2i, 3)), ABS(UBOUND(c2i, 3)))
     484           10 :       ALLOCATE (cell_to_index(-nx:nx, -ny:ny, -nz:nz))
     485              :       cell_to_index(LBOUND(c2i, 1):UBOUND(c2i, 1), &
     486              :                     LBOUND(c2i, 2):UBOUND(c2i, 2), &
     487           84 :                     LBOUND(c2i, 3):UBOUND(c2i, 3)) = c2i
     488              : 
     489              :       ! identify cells with no mirror img
     490              :       nomirror = 0
     491           52 :       DO ic = 1, ncell
     492          200 :          cell = i2c(:, ic)
     493           50 :          IF (cell_to_index(-cell(1), -cell(2), -cell(3)) == 0) &
     494            6 :             nomirror = nomirror + 1
     495              :       END DO
     496              : 
     497              :       ! create the mirror imgs
     498            6 :       ALLOCATE (index_to_cell(3, ncell + nomirror))
     499          202 :       index_to_cell(:, 1:ncell) = i2c
     500              : 
     501              :       nomirror = 0 ! count the imgs without mirror
     502           52 :       DO ic = 1, ncell
     503          200 :          cell = index_to_cell(:, ic)
     504           52 :          IF (cell_to_index(-cell(1), -cell(2), -cell(3)) == 0) THEN
     505            4 :             nomirror = nomirror + 1
     506           16 :             index_to_cell(:, ncell + nomirror) = -cell
     507            4 :             cell_to_index(-cell(1), -cell(2), -cell(3)) = ncell + nomirror
     508              :          END IF
     509              :       END DO
     510            2 :       ncell = ncell + nomirror
     511              : 
     512            2 :       CALL get_kpoint_info(kpoints, sab_nl=sab_nl)
     513              :       ! allocate the nonsymmetric matrices
     514            2 :       NULLIFY (mat_nosym)
     515            2 :       CALL dbcsr_allocate_matrix_set(mat_nosym, nspin, ncell)
     516            4 :       DO ispin = 1, nspin
     517           58 :          DO ic = 1, ncell
     518           54 :             ALLOCATE (mat_nosym(ispin, ic)%matrix)
     519              :             CALL dbcsr_create(matrix=mat_nosym(ispin, ic)%matrix, &
     520              :                               template=mat(1, 1)%matrix, &
     521           54 :                               matrix_type=dbcsr_type_no_symmetry)
     522              :             CALL cp_dbcsr_alloc_block_from_nbl(mat_nosym(ispin, ic)%matrix, &
     523           54 :                                                sab_nl, desymmetrize=.TRUE.)
     524           56 :             CALL dbcsr_set(mat_nosym(ispin, ic)%matrix, 0.0_dp)
     525              :          END DO
     526              :       END DO
     527              : 
     528            4 :       DO ispin = 1, nspin
     529              :          ! desymmetrize the matrix for real space printing
     530            2 :          CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
     531          506 :          DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
     532          504 :             CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom, cell=cell)
     533              : 
     534          504 :             ic = cell_to_index(cell(1), cell(2), cell(3))
     535          504 :             icn = cell_to_index(-cell(1), -cell(2), -cell(3))
     536          504 :             CPASSERT(icn > 0)
     537              : 
     538          504 :             irow = iatom
     539          504 :             icol = jatom
     540          504 :             lwtr = .FALSE.
     541              :             ! always copy from the top
     542          504 :             IF (iatom > jatom) THEN
     543          216 :                irow = jatom
     544          216 :                icol = iatom
     545          216 :                lwtr = .TRUE.
     546              :             END IF
     547              : 
     548              :             CALL dbcsr_get_block_p(matrix=mat(ispin, ic)%matrix, &
     549          504 :                                    row=irow, col=icol, block=block, found=found)
     550          504 :             CPASSERT(found)
     551              : 
     552              :             ! copy to M(R) at (iatom,jatom)
     553              :             ! copy to M(-R) at (jatom,iatom)
     554          506 :             IF (lwtr) THEN
     555              :                CALL dbcsr_put_block(matrix=mat_nosym(ispin, ic)%matrix, &
     556         4536 :                                     row=iatom, col=jatom, block=TRANSPOSE(block))
     557              :                CALL dbcsr_put_block(matrix=mat_nosym(ispin, icn)%matrix, &
     558          216 :                                     row=jatom, col=iatom, block=block)
     559              :             ELSE
     560              :                CALL dbcsr_put_block(matrix=mat_nosym(ispin, ic)%matrix, &
     561          288 :                                     row=iatom, col=jatom, block=block)
     562              :                CALL dbcsr_put_block(matrix=mat_nosym(ispin, icn)%matrix, &
     563         6048 :                                     row=jatom, col=iatom, block=TRANSPOSE(block))
     564              :             END IF
     565              :          END DO
     566            4 :          CALL neighbor_list_iterator_release(nl_iterator)
     567              :       END DO
     568              : 
     569            4 :       DO ispin = 1, nspin
     570           58 :          DO ic = 1, ncell
     571           56 :             CALL dbcsr_finalize(mat_nosym(ispin, ic)%matrix)
     572              :          END DO
     573              :       END DO
     574              : 
     575            2 :       CALL timestop(handle)
     576              : 
     577            2 :    END SUBROUTINE desymmetrize_rs_matrix
     578              : 
     579              : END MODULE qs_scf_csr_write
        

Generated by: LCOV version 2.0-1