comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-28 22:25:32 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Enable printing of the real space representation of the KS and S matrices in runs with kpoints

Browse source
This commit is contained in:
Fabian Ducry 2021-04-30 15:05:07 +02:00 • committed by GitHub
parent 63dfad33bf
commit d8b6a4f63f
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
5 changed files with 483 additions and 141 deletions
Download patch file Download diff file Expand all files Collapse all files

16
src/input_cp2k_dft.F
View file

@ -2255,6 +2255,14 @@ CONTAINS
lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="Real_space", &
description="Print the KS matrix in real-space instead of k-space.. ", &
usage="REAL_SPACE {logical}", &
repeats=.FALSE., &
default_l_val=.FALSE., &
lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL section_add_subsection(section, print_key)
CALL section_release(print_key)
@ -2287,6 +2295,14 @@ CONTAINS
lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="Real_space", &
description="Print the overlap matrix in real-space instead of k-space.. ", &
usage="REAL_SPACE {logical}", &
repeats=.FALSE., &
default_l_val=.FALSE., &
lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL section_add_subsection(section, print_key)
CALL section_release(print_key)

542
src/qs_scf_csr_write.F
View file

@ -13,6 +13,8 @@
! **************************************************************************************************
MODULE qs_scf_csr_write
USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,&
dbcsr_deallocate_matrix_set
USE cp_log_handling, ONLY: cp_get_default_logger,&
cp_logger_get_default_io_unit,&
cp_logger_type
@ -23,9 +25,10 @@ MODULE qs_scf_csr_write
USE dbcsr_api, ONLY: &
dbcsr_add, dbcsr_convert_dbcsr_to_csr, dbcsr_copy, dbcsr_create, &
dbcsr_csr_create_from_dbcsr, dbcsr_csr_dbcsr_blkrow_dist, dbcsr_csr_destroy, &
dbcsr_csr_type, dbcsr_csr_write, dbcsr_desymmetrize, dbcsr_has_symmetry, dbcsr_p_type, &
dbcsr_release, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_complex_8, &
dbcsr_type_no_symmetry, dbcsr_type_symmetric
dbcsr_csr_type, dbcsr_csr_write, dbcsr_desymmetrize, dbcsr_finalize, dbcsr_get_block_p, &
dbcsr_has_symmetry, dbcsr_p_type, dbcsr_put_block, dbcsr_release, dbcsr_set, dbcsr_type, &
dbcsr_type_antisymmetric, dbcsr_type_complex_8, dbcsr_type_no_symmetry, dbcsr_type_real_8, &
dbcsr_type_symmetric
USE input_section_types, ONLY: section_vals_get_subs_vals,&
section_vals_type,&
section_vals_val_get
@ -36,7 +39,13 @@ MODULE qs_scf_csr_write
kpoint_type
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
USE qs_neighbor_list_types, ONLY: get_iterator_info,&
get_neighbor_list_set_p,&
neighbor_list_iterate,&
neighbor_list_iterator_create,&
neighbor_list_iterator_p_type,&
neighbor_list_iterator_release,&
neighbor_list_set_p_type
#include "./base/base_uses.f90"
IMPLICIT NONE
@ -66,10 +75,9 @@ CONTAINS
CHARACTER(len=*), PARAMETER :: routineN = 'write_ks_matrix_csr'
INTEGER :: handle, output_unit
LOGICAL :: do_kpoints, do_ks_csr_write
LOGICAL :: do_kpoints, do_ks_csr_write, real_space
TYPE(cp_logger_type), POINTER :: logger
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks
TYPE(kpoint_type), POINTER :: kpoints
TYPE(section_vals_type), POINTER :: dft_section
@ -85,15 +93,15 @@ CONTAINS
"PRINT%KS_CSR_WRITE"), cp_p_file)
IF (do_ks_csr_write) THEN
CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints)
CALL get_qs_env(qs_env=qs_env, kpoints=kpoints, matrix_ks_kp=matrix_ks, do_kpoints=do_kpoints)
CALL section_vals_val_get(dft_section, "PRINT%KS_CSR_WRITE%REAL_SPACE", &
l_val=real_space)
IF (do_kpoints) THEN
CALL get_qs_env(qs_env=qs_env, kpoints=kpoints, matrix_ks_kp=matrix_ks_kp)
CALL write_matrix_csr(dft_section, matrix_to_print_kp=matrix_ks_kp, &
kpoints=kpoints, prefix="KS")
IF (do_kpoints .AND. .NOT. real_space) THEN
CALL write_matrix_kp_csr(mat=matrix_ks, dft_section=dft_section, &
kpoints=kpoints, prefix="KS")
ELSE
CALL get_qs_env(qs_env=qs_env, matrix_ks=matrix_ks)
CALL write_matrix_csr(dft_section, matrix_to_print=matrix_ks, &
CALL write_matrix_csr(dft_section, mat=matrix_ks, kpoints=kpoints, do_kpoints=do_kpoints, &
prefix="KS")
END IF
END IF
@ -117,10 +125,9 @@ CONTAINS
CHARACTER(len=*), PARAMETER :: routineN = 'write_s_matrix_csr'
INTEGER :: handle, output_unit
LOGICAL :: do_kpoints, do_s_csr_write
LOGICAL :: do_kpoints, do_s_csr_write, real_space
TYPE(cp_logger_type), POINTER :: logger
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_kp
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
TYPE(kpoint_type), POINTER :: kpoints
TYPE(section_vals_type), POINTER :: dft_section
@ -136,15 +143,15 @@ CONTAINS
"PRINT%S_CSR_WRITE"), cp_p_file)
IF (do_s_csr_write) THEN
CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints)
CALL get_qs_env(qs_env=qs_env, kpoints=kpoints, matrix_s_kp=matrix_s, do_kpoints=do_kpoints)
CALL section_vals_val_get(dft_section, "PRINT%S_CSR_WRITE%REAL_SPACE", &
l_val=real_space)
IF (do_kpoints) THEN
CALL get_qs_env(qs_env=qs_env, kpoints=kpoints, matrix_s_kp=matrix_s_kp)
CALL write_matrix_csr(dft_section, matrix_to_print_kp=matrix_s_kp, &
kpoints=kpoints, prefix="S")
IF (do_kpoints .AND. .NOT. real_space) THEN
CALL write_matrix_kp_csr(mat=matrix_s, dft_section=dft_section, &
kpoints=kpoints, prefix="S")
ELSE
CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
CALL write_matrix_csr(dft_section, matrix_to_print=matrix_s, &
CALL write_matrix_csr(dft_section, mat=matrix_s, kpoints=kpoints, do_kpoints=do_kpoints, &
prefix="S")
END IF
END IF
@ -154,25 +161,165 @@ CONTAINS
END SUBROUTINE write_s_matrix_csr
! **************************************************************************************************
!> \brief helper function to print the KS or S matrix to file
!> \brief helper function to print the real space representation of KS or S matrix to file
!> \param dft_section the dft_section
!> \param matrix_to_print Hamiltonian or overlap matrix for Gamma point calculations
!> \param matrix_to_print_kp Hamiltonian or overlap matrix for k-point calculations
!> \param kpoints the kpoints
!> \param mat Hamiltonian or overlap matrix
!> \param kpoints Kpoint environment
!> \param prefix string to distinguish between KS and S matrix
!> \param do_kpoints Whether it is a gamma-point or k-point calculation
!> \par History
!> Moved most of the code from write_ks_matrix_csr and write_s_matrix_csr
!> Removed the code for printing k-point dependent matrices and added
!> printing of the real space representation
! **************************************************************************************************
SUBROUTINE write_matrix_csr(dft_section, matrix_to_print, matrix_to_print_kp, kpoints, prefix)
SUBROUTINE write_matrix_csr(dft_section, mat, kpoints, prefix, do_kpoints)
TYPE(section_vals_type), INTENT(IN), POINTER :: dft_section
TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN), &
OPTIONAL, POINTER :: matrix_to_print
TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(IN), &
OPTIONAL, POINTER :: matrix_to_print_kp
TYPE(kpoint_type), INTENT(IN), OPTIONAL, POINTER :: kpoints
POINTER :: mat
TYPE(kpoint_type), INTENT(IN), POINTER :: kpoints
CHARACTER(*), INTENT(in) :: prefix
LOGICAL, INTENT(IN) :: do_kpoints
CHARACTER(len=*), PARAMETER :: routineN = 'write_matrix_csr'
CHARACTER(LEN=default_path_length) :: file_name, fileformat, subs_string
INTEGER :: handle, ic, ispin, ncell, nspin, &
output_unit, unit_nr
INTEGER, ALLOCATABLE, DIMENSION(:, :) :: index_to_cell
INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cell_to_index
INTEGER, DIMENSION(:, :), POINTER :: i2c
INTEGER, DIMENSION(:, :, :), POINTER :: c2i
LOGICAL :: bin, do_symmetric, uptr
REAL(KIND=dp) :: thld
TYPE(cp_logger_type), POINTER :: logger
TYPE(dbcsr_csr_type) :: mat_csr
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_nosym
TYPE(dbcsr_type) :: matrix_nosym
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_nl
CALL timeset(routineN, handle)
logger => cp_get_default_logger()
output_unit = cp_logger_get_default_io_unit(logger)
subs_string = "PRINT%"//prefix//"_CSR_WRITE"
CALL section_vals_val_get(dft_section, subs_string//"%THRESHOLD", r_val=thld)
CALL section_vals_val_get(dft_section, subs_string//"%UPPER_TRIANGULAR", l_val=uptr)
CALL section_vals_val_get(dft_section, subs_string//"%BINARY", l_val=bin)
IF (bin) THEN
fileformat = "UNFORMATTED"
ELSE
fileformat = "FORMATTED"
END IF
nspin = SIZE(mat, 1)
ncell = SIZE(mat, 2)
IF (do_kpoints) THEN
i2c => kpoints%index_to_cell
c2i => kpoints%cell_to_index
NULLIFY (sab_nl)
CALL get_kpoint_info(kpoints, sab_nl=sab_nl)
CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
! desymmetrize the KS or S matrices if necessary
IF (do_symmetric) THEN
CALL desymmetrize_rs_matrix(mat, mat_nosym, cell_to_index, index_to_cell, kpoints)
ncell = SIZE(index_to_cell, 2) ! update the number of cells
ELSE
ALLOCATE (cell_to_index(LBOUND(c2i, 1):UBOUND(c2i, 1), &
LBOUND(c2i, 2):UBOUND(c2i, 2), &
LBOUND(c2i, 3):UBOUND(c2i, 3)))
cell_to_index(LBOUND(c2i, 1):UBOUND(c2i, 1), &
LBOUND(c2i, 2):UBOUND(c2i, 2), &
LBOUND(c2i, 3):UBOUND(c2i, 3)) = c2i
ALLOCATE (index_to_cell(3, ncell))
index_to_cell(1:3, 1:ncell) = i2c
mat_nosym => mat
END IF
! print the index to cell mapping to the output
IF (output_unit > 0) THEN
WRITE (output_unit, "(/,A15,T15,I4,A)") prefix//" CSR write| ", &
ncell, " periodic images"
WRITE (output_unit, "(T7,A,T17,A,T24,A,T31,A)") "Number", "X", "Y", "Z"
DO ic = 1, ncell
WRITE (output_unit, "(T8,I3,T15,I3,T22,I3,T29,I3)") ic, index_to_cell(:, ic)
END DO
END IF
END IF
! write the csr file(s)
DO ispin = 1, nspin
DO ic = 1, ncell
IF (do_kpoints) THEN
CALL dbcsr_copy(matrix_nosym, mat_nosym(ispin, ic)%matrix)
WRITE (file_name, '(2(A,I0))') prefix//"_SPIN_", ispin, "_R_", ic
ELSE
IF (dbcsr_has_symmetry(mat(ispin, ic)%matrix)) THEN
CALL dbcsr_desymmetrize(mat(ispin, ic)%matrix, matrix_nosym)
ELSE
CALL dbcsr_copy(matrix_nosym, mat(ispin, ic)%matrix)
END IF
WRITE (file_name, '(A,I0)') prefix//"_SPIN_", ispin
END IF
! Convert dbcsr to csr
CALL dbcsr_csr_create_from_dbcsr(matrix_nosym, &
mat_csr, dbcsr_csr_dbcsr_blkrow_dist)
CALL dbcsr_convert_dbcsr_to_csr(matrix_nosym, mat_csr)
! Write to file
unit_nr = cp_print_key_unit_nr(logger, dft_section, subs_string, &
extension=".csr", middle_name=TRIM(file_name), &
file_status="REPLACE", file_form=fileformat)
CALL dbcsr_csr_write(mat_csr, unit_nr, upper_triangle=uptr, threshold=thld, binary=bin)
CALL cp_print_key_finished_output(unit_nr, logger, dft_section, subs_string)
CALL dbcsr_csr_destroy(mat_csr)
CALL dbcsr_release(matrix_nosym)
END DO
END DO
! clean up
IF (do_kpoints) THEN
DEALLOCATE (cell_to_index, index_to_cell)
IF (do_symmetric) THEN
DO ispin = 1, nspin
DO ic = 1, ncell
CALL dbcsr_release(mat_nosym(ispin, ic)%matrix)
END DO
END DO
CALL dbcsr_deallocate_matrix_set(mat_nosym)
END IF
END IF
CALL timestop(handle)
END SUBROUTINE write_matrix_csr
! **************************************************************************************************
!> \brief helper function to print the k-dependent KS or S matrix to file
!> \param mat Hamiltonian or overlap matrix for k-point calculations
!> \param dft_section the dft_section
!> \param kpoints Kpoint environment
!> \param prefix string to distinguish between KS and S matrix
!> \par History
!> Moved the code from write_matrix_csr to write_matrix_kp_csr
!> \author Fabian Ducry
! **************************************************************************************************
SUBROUTINE write_matrix_kp_csr(mat, dft_section, kpoints, prefix)
TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(IN), &
POINTER :: mat
TYPE(section_vals_type), INTENT(IN), POINTER :: dft_section
TYPE(kpoint_type), INTENT(IN), POINTER :: kpoints
CHARACTER(*), INTENT(in) :: prefix
CHARACTER(len=*), PARAMETER :: routineN = 'write_matrix_kp_csr'
COMPLEX(KIND=dp), PARAMETER :: cone = CMPLX(1.0_dp, 0.0_dp, KIND=dp), &
ione = CMPLX(0.0_dp, 1.0_dp, KIND=dp)
@ -181,7 +328,6 @@ CONTAINS
nkp_groups, output_unit, unit_nr
INTEGER, DIMENSION(2) :: kp_range
INTEGER, DIMENSION(:, :), POINTER :: kp_dist
INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
LOGICAL :: bin, uptr, use_real_wfn
REAL(KIND=dp) :: thld
REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
@ -210,122 +356,238 @@ CONTAINS
fileformat = "FORMATTED"
END IF
IF (PRESENT(kpoints)) THEN
! Calculate the Hamiltonian at the k-points
NULLIFY (sab_nl)
CALL get_kpoint_info(kpoints, xkp=xkp, use_real_wfn=use_real_wfn, kp_range=kp_range, &
nkp_groups=nkp_groups, kp_dist=kp_dist, sab_nl=sab_nl, &
cell_to_index=cell_to_index)
NULLIFY (sab_nl)
ALLOCATE (rmatrix)
CALL dbcsr_create(rmatrix, template=matrix_to_print_kp(1, 1)%matrix, &
matrix_type=dbcsr_type_symmetric)
CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_nl)
! Calculate the Hamiltonian at the k-points
CALL get_kpoint_info(kpoints, xkp=xkp, use_real_wfn=use_real_wfn, kp_range=kp_range, &
nkp_groups=nkp_groups, kp_dist=kp_dist, sab_nl=sab_nl)
IF (.NOT. use_real_wfn) THEN
! Allocate temporary variables
ALLOCATE (rmatrix_nosym, imatrix, imatrix_nosym, cmatrix, tmatrix)
CALL dbcsr_create(rmatrix_nosym, template=matrix_to_print_kp(1, 1)%matrix, &
matrix_type=dbcsr_type_no_symmetry)
CALL dbcsr_create(imatrix, template=matrix_to_print_kp(1, 1)%matrix, &
matrix_type=dbcsr_type_antisymmetric)
CALL dbcsr_create(imatrix_nosym, template=matrix_to_print_kp(1, 1)%matrix, &
matrix_type=dbcsr_type_no_symmetry)
CALL dbcsr_create(cmatrix, template=matrix_to_print_kp(1, 1)%matrix, &
matrix_type=dbcsr_type_no_symmetry, &
data_type=dbcsr_type_complex_8)
CALL dbcsr_create(tmatrix, template=matrix_to_print_kp(1, 1)%matrix, &
matrix_type=dbcsr_type_no_symmetry, &
data_type=dbcsr_type_complex_8)
CALL cp_dbcsr_alloc_block_from_nbl(rmatrix_nosym, sab_nl)
CALL cp_dbcsr_alloc_block_from_nbl(imatrix, sab_nl)
CALL cp_dbcsr_alloc_block_from_nbl(imatrix_nosym, sab_nl)
CALL cp_dbcsr_alloc_block_from_nbl(cmatrix, sab_nl)
CALL cp_dbcsr_alloc_block_from_nbl(tmatrix, sab_nl)
END IF
ALLOCATE (rmatrix)
CALL dbcsr_create(rmatrix, template=mat(1, 1)%matrix, &
matrix_type=dbcsr_type_symmetric)
CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_nl)
kplocal = kp_range(2) - kp_range(1) + 1
DO ikp = 1, kplocal
DO ispin = 1, SIZE(matrix_to_print_kp, 1)
DO igroup = 1, nkp_groups
! number of current kpoint
ik = kp_dist(1, igroup) + ikp - 1
CALL dbcsr_set(rmatrix, 0.0_dp)
IF (use_real_wfn) THEN
! FT of KS matrix
CALL rskp_transform(rmatrix=rmatrix, rsmat=matrix_to_print_kp, ispin=ispin, &
xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
! Convert to desymmetrized csr matrix
CALL dbcsr_desymmetrize(rmatrix, matrix_nosym)
CALL dbcsr_csr_create_from_dbcsr(matrix_nosym, mat_csr, dbcsr_csr_dbcsr_blkrow_dist)
CALL dbcsr_convert_dbcsr_to_csr(matrix_nosym, mat_csr)
CALL dbcsr_release(matrix_nosym)
ELSE
! FT of KS matrix
CALL dbcsr_set(imatrix, 0.0_dp)
CALL rskp_transform(rmatrix=rmatrix, cmatrix=imatrix, rsmat=matrix_to_print_kp, ispin=ispin, &
xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
IF (.NOT. use_real_wfn) THEN
! Allocate temporary variables
ALLOCATE (rmatrix_nosym, imatrix, imatrix_nosym, cmatrix, tmatrix)
CALL dbcsr_create(rmatrix_nosym, template=mat(1, 1)%matrix, &
matrix_type=dbcsr_type_no_symmetry)
CALL dbcsr_create(imatrix, template=mat(1, 1)%matrix, &
matrix_type=dbcsr_type_antisymmetric)
CALL dbcsr_create(imatrix_nosym, template=mat(1, 1)%matrix, &
matrix_type=dbcsr_type_no_symmetry)
CALL dbcsr_create(cmatrix, template=mat(1, 1)%matrix, &
matrix_type=dbcsr_type_no_symmetry, &
data_type=dbcsr_type_complex_8)
CALL dbcsr_create(tmatrix, template=mat(1, 1)%matrix, &
matrix_type=dbcsr_type_no_symmetry, &
data_type=dbcsr_type_complex_8)
CALL cp_dbcsr_alloc_block_from_nbl(rmatrix_nosym, sab_nl)
CALL cp_dbcsr_alloc_block_from_nbl(imatrix, sab_nl)
CALL cp_dbcsr_alloc_block_from_nbl(imatrix_nosym, sab_nl)
CALL cp_dbcsr_alloc_block_from_nbl(cmatrix, sab_nl)
CALL cp_dbcsr_alloc_block_from_nbl(tmatrix, sab_nl)
END IF
! Desymmetrize and sum the real and imaginary part into
! cmatrix
CALL dbcsr_desymmetrize(rmatrix, rmatrix_nosym)
CALL dbcsr_desymmetrize(imatrix, imatrix_nosym)
CALL dbcsr_copy(cmatrix, rmatrix_nosym)
CALL dbcsr_copy(tmatrix, imatrix_nosym)
CALL dbcsr_add(cmatrix, tmatrix, cone, ione)
! Convert to csr
CALL dbcsr_csr_create_from_dbcsr(cmatrix, mat_csr, dbcsr_csr_dbcsr_blkrow_dist)
CALL dbcsr_convert_dbcsr_to_csr(cmatrix, mat_csr)
END IF
! Write to file
WRITE (file_name, '(2(A,I0))') prefix//"_SPIN_", ispin, "_K_", ik
unit_nr = cp_print_key_unit_nr(logger, dft_section, subs_string, &
extension=".csr", middle_name=TRIM(file_name), &
file_status="REPLACE", file_form=fileformat)
CALL dbcsr_csr_write(mat_csr, unit_nr, upper_triangle=uptr, threshold=thld, binary=bin)
kplocal = kp_range(2) - kp_range(1) + 1
DO ikp = 1, kplocal
DO ispin = 1, SIZE(mat, 1)
DO igroup = 1, nkp_groups
! number of current kpoint
ik = kp_dist(1, igroup) + ikp - 1
CALL dbcsr_set(rmatrix, 0.0_dp)
IF (use_real_wfn) THEN
! FT of the matrix
CALL rskp_transform(rmatrix=rmatrix, rsmat=mat, ispin=ispin, &
xkp=xkp(1:3, ik), cell_to_index=kpoints%cell_to_index, sab_nl=sab_nl)
! Convert to desymmetrized csr matrix
CALL dbcsr_desymmetrize(rmatrix, matrix_nosym)
CALL dbcsr_csr_create_from_dbcsr(matrix_nosym, mat_csr, dbcsr_csr_dbcsr_blkrow_dist)
CALL dbcsr_convert_dbcsr_to_csr(matrix_nosym, mat_csr)
CALL dbcsr_release(matrix_nosym)
ELSE
! FT of the matrix
CALL dbcsr_set(imatrix, 0.0_dp)
CALL rskp_transform(rmatrix=rmatrix, cmatrix=imatrix, rsmat=mat, ispin=ispin, &
xkp=xkp(1:3, ik), cell_to_index=kpoints%cell_to_index, sab_nl=sab_nl)
CALL cp_print_key_finished_output(unit_nr, logger, dft_section, subs_string)
! Desymmetrize and sum the real and imaginary part into
! cmatrix
CALL dbcsr_desymmetrize(rmatrix, rmatrix_nosym)
CALL dbcsr_desymmetrize(imatrix, imatrix_nosym)
CALL dbcsr_copy(cmatrix, rmatrix_nosym)
CALL dbcsr_copy(tmatrix, imatrix_nosym)
CALL dbcsr_add(cmatrix, tmatrix, cone, ione)
! Convert to csr
CALL dbcsr_csr_create_from_dbcsr(cmatrix, mat_csr, dbcsr_csr_dbcsr_blkrow_dist)
CALL dbcsr_convert_dbcsr_to_csr(cmatrix, mat_csr)
END IF
! Write to file
WRITE (file_name, '(2(A,I0))') prefix//"_SPIN_", ispin, "_K_", ik
unit_nr = cp_print_key_unit_nr(logger, dft_section, subs_string, &
extension=".csr", middle_name=TRIM(file_name), &
file_status="REPLACE", file_form=fileformat)
CALL dbcsr_csr_write(mat_csr, unit_nr, upper_triangle=uptr, threshold=thld, binary=bin)
CALL dbcsr_csr_destroy(mat_csr)
END DO
CALL cp_print_key_finished_output(unit_nr, logger, dft_section, subs_string)
CALL dbcsr_csr_destroy(mat_csr)
END DO
END DO
CALL dbcsr_release(rmatrix)
DEALLOCATE (rmatrix)
IF (.NOT. use_real_wfn) THEN
CALL dbcsr_release(rmatrix_nosym)
CALL dbcsr_release(imatrix)
CALL dbcsr_release(imatrix_nosym)
CALL dbcsr_release(cmatrix)
CALL dbcsr_release(tmatrix)
DEALLOCATE (rmatrix_nosym, imatrix, imatrix_nosym, cmatrix, tmatrix)
END IF
ELSE
DO ispin = 1, SIZE(matrix_to_print)
! Desymmetrize
IF (dbcsr_has_symmetry(matrix_to_print(ispin)%matrix)) THEN
CALL dbcsr_desymmetrize(matrix_to_print(ispin)%matrix, matrix_nosym)
ELSE
CALL dbcsr_copy(matrix_nosym, matrix_to_print(ispin)%matrix)
END IF
! Convert dbcsr to csr
CALL dbcsr_csr_create_from_dbcsr(matrix_nosym, mat_csr, dbcsr_csr_dbcsr_blkrow_dist)
CALL dbcsr_convert_dbcsr_to_csr(matrix_nosym, mat_csr)
! Write to file
WRITE (file_name, '(A,I0)') prefix//"_SPIN_", ispin
unit_nr = cp_print_key_unit_nr(logger, dft_section, subs_string, &
extension=".csr", middle_name=TRIM(file_name), &
file_status="REPLACE", file_form=fileformat)
CALL dbcsr_csr_write(mat_csr, unit_nr, upper_triangle=uptr, threshold=thld, binary=bin)
CALL cp_print_key_finished_output(unit_nr, logger, dft_section, subs_string)
CALL dbcsr_csr_destroy(mat_csr)
CALL dbcsr_release(matrix_nosym)
END DO
END DO
CALL dbcsr_release(rmatrix)
DEALLOCATE (rmatrix)
IF (.NOT. use_real_wfn) THEN
CALL dbcsr_release(rmatrix_nosym)
CALL dbcsr_release(imatrix)
CALL dbcsr_release(imatrix_nosym)
CALL dbcsr_release(cmatrix)
CALL dbcsr_release(tmatrix)
DEALLOCATE (rmatrix_nosym, imatrix, imatrix_nosym, cmatrix, tmatrix)
END IF
CALL timestop(handle)
END SUBROUTINE write_matrix_csr
END SUBROUTINE write_matrix_kp_csr
! **************************************************************************************************
!> \brief Desymmetrizes the KS or S matrices which are stored in symmetric !matrices
!> \param mat Hamiltonian or overlap matrices
!> \param mat_nosym Desymmetrized Hamiltonian or overlap matrices
!> \param cell_to_index Mapping of cell indices to linear RS indices
!> \param index_to_cell Mapping of linear RS indices to cell indices
!> \param kpoints Kpoint environment
!> \author Fabian Ducry
! **************************************************************************************************
SUBROUTINE desymmetrize_rs_matrix(mat, mat_nosym, cell_to_index, index_to_cell, kpoints)
TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(IN), &
POINTER :: mat
TYPE(dbcsr_p_type), DIMENSION(:, :), &
INTENT(INOUT), POINTER :: mat_nosym
INTEGER, ALLOCATABLE, DIMENSION(:, :, :), &
INTENT(OUT) :: cell_to_index
INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(OUT) :: index_to_cell
TYPE(kpoint_type), INTENT(IN), POINTER :: kpoints
CHARACTER(len=*), PARAMETER :: routineN = 'desymmetrize_rs_matrix'
INTEGER :: handle, iatom, ic, icn, icol, irow, &
ispin, jatom, ncell, nomirror, nspin, &
nx, ny, nz
INTEGER, DIMENSION(3) :: cell
INTEGER, DIMENSION(:, :), POINTER :: i2c
INTEGER, DIMENSION(:, :, :), POINTER :: c2i
LOGICAL :: found, lwtr
REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
TYPE(neighbor_list_iterator_p_type), &
DIMENSION(:), POINTER :: nl_iterator
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_nl
CALL timeset(routineN, handle)
i2c => kpoints%index_to_cell
c2i => kpoints%cell_to_index
ncell = SIZE(i2c, 2)
nspin = SIZE(mat, 1)
nx = MAX(ABS(LBOUND(c2i, 1)), ABS(UBOUND(c2i, 1)))
ny = MAX(ABS(LBOUND(c2i, 2)), ABS(UBOUND(c2i, 3)))
nz = MAX(ABS(LBOUND(c2i, 3)), ABS(UBOUND(c2i, 3)))
ALLOCATE (cell_to_index(-nx:nx, -ny:ny, -nz:nz))
cell_to_index(LBOUND(c2i, 1):UBOUND(c2i, 1), &
LBOUND(c2i, 2):UBOUND(c2i, 2), &
LBOUND(c2i, 3):UBOUND(c2i, 3)) = c2i
! identify cells with no mirror img
nomirror = 0
DO ic = 1, ncell
cell = i2c(:, ic)
IF (cell_to_index(-cell(1), -cell(2), -cell(3)) == 0) &
nomirror = nomirror + 1
END DO
! create the mirror imgs
ALLOCATE (index_to_cell(3, ncell + nomirror))
index_to_cell(:, 1:ncell) = i2c
nomirror = 0 ! count the imgs without mirror
DO ic = 1, ncell
cell = index_to_cell(:, ic)
IF (cell_to_index(-cell(1), -cell(2), -cell(3)) == 0) THEN
nomirror = nomirror + 1
index_to_cell(:, ncell + nomirror) = -cell
cell_to_index(-cell(1), -cell(2), -cell(3)) = ncell + nomirror
END IF
END DO
ncell = ncell + nomirror
CALL get_kpoint_info(kpoints, sab_nl=sab_nl)
! allocate the nonsymmetric matrices
NULLIFY (mat_nosym)
CALL dbcsr_allocate_matrix_set(mat_nosym, nspin, ncell)
DO ispin = 1, nspin
DO ic = 1, ncell
ALLOCATE (mat_nosym(ispin, ic)%matrix)
CALL dbcsr_create(matrix=mat_nosym(ispin, ic)%matrix, &
template=mat(1, 1)%matrix, &
matrix_type=dbcsr_type_no_symmetry, &
data_type=dbcsr_type_real_8)
CALL cp_dbcsr_alloc_block_from_nbl(mat_nosym(ispin, ic)%matrix, &
sab_nl, desymmetrize=.TRUE.)
CALL dbcsr_set(mat_nosym(ispin, ic)%matrix, 0.0_dp)
END DO
END DO
DO ispin = 1, nspin
! desymmetrize the matrix for real space printing
CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom, cell=cell)
ic = cell_to_index(cell(1), cell(2), cell(3))
icn = cell_to_index(-cell(1), -cell(2), -cell(3))
CPASSERT(icn > 0)
irow = iatom
icol = jatom
lwtr = .FALSE.
! always copy from the top
IF (iatom > jatom) THEN
irow = jatom
icol = iatom
lwtr = .TRUE.
END IF
CALL dbcsr_get_block_p(matrix=mat(ispin, ic)%matrix, &
row=irow, col=icol, block=block, found=found)
CPASSERT(found)
! copy to M(R) at (iatom,jatom)
! copy to M(-R) at (jatom,iatom)
IF (lwtr) THEN
CALL dbcsr_put_block(matrix=mat_nosym(ispin, ic)%matrix, &
row=iatom, col=jatom, block=TRANSPOSE(block))
CALL dbcsr_put_block(matrix=mat_nosym(ispin, icn)%matrix, &
row=jatom, col=iatom, block=block)
ELSE
CALL dbcsr_put_block(matrix=mat_nosym(ispin, ic)%matrix, &
row=iatom, col=jatom, block=block)
CALL dbcsr_put_block(matrix=mat_nosym(ispin, icn)%matrix, &
row=jatom, col=iatom, block=TRANSPOSE(block))
END IF
END DO
CALL neighbor_list_iterator_release(nl_iterator)
END DO
DO ispin = 1, nspin
DO ic = 1, ncell
CALL dbcsr_finalize(mat_nosym(ispin, ic)%matrix)
END DO
END DO
CALL timestop(handle)
END SUBROUTINE desymmetrize_rs_matrix
END MODULE qs_scf_csr_write

1
tests/QS/regtest-kp-3/TEST_FILES
View file

@ -4,6 +4,7 @@
# 1 compares the last total energy in the file
# for details see cp2k/tools/do_regtest
cc1.inp 1 1.0E-12 -45.38964266770495
bn1.inp 94 1.0E-14 27
#cc2.inp 1 1.0E-12 -45.38989062944852
#cc3.inp 1 1.0E-12 -45.38990176407365
#cc4.inp 1 1.0E-12 -45.38975726647867

62
tests/QS/regtest-kp-3/bn1.inp Normal file
View file

@ -0,0 +1,62 @@
&FORCE_EVAL
&DFT
&QS
METHOD GPW
&END QS
&KPOINTS
SCHEME MONKHORST-PACK 2 2 1
&END KPOINTS
&XC
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&SCF
SCF_GUESS ATOMIC
MAX_SCF 1
&MIXING
METHOD DIRECT_P_MIXING
ALPHA 0.40
&END
&END SCF
&PRINT
&KS_CSR_WRITE
BINARY .TRUE.
THRESHOLD 1.0E-4
REAL_SPACE .TRUE.
&END
&END
BASIS_SET_FILE_NAME BASIS_MOLOPT
&END DFT
&SUBSYS
&KIND B
BASIS_SET SZV-MOLOPT-SR-GTH
POTENTIAL GTH-PBE
&END KIND
&KIND N
BASIS_SET SZV-MOLOPT-SR-GTH
POTENTIAL GTH-PBE
&END KIND
&TOPOLOGY
&END
&CELL
ABC 4.3483 5.0209 15.0
&END CELL
&COORD
SCALED
N 0.0000 0 0
B 0.3333 0 0
N 0.5000 0.2500 0
B 0.8333 0.2500 0
N 0.0000 0.5000 0
B 0.3333 0.5000 0
N 0.5000 0.7500 0
B 0.8333 0.7500 0
&END COORD
&END SUBSYS
STRESS_TENSOR ANALYTICAL
&END FORCE_EVAL
&GLOBAL
PROJECT BN
RUN_TYPE energy
PRINT_LEVEL MEDIUM
&END GLOBAL

3
tests/TEST_TYPES
View file

@ -1,4 +1,4 @@
93
94
Total energy:!3
MD| Potential energy!5
Total energy \[eV\]:!4
@ -92,6 +92,7 @@ Final localization: !3
Ionization potentials for XPS !8
FCIDUMP| Checksum: !3
SPGR| SPACE GROUP NUMBER: !5
KS CSR write| !4
#
# these are the tests the can be selected for regtesting.
# do regtest will grep for test_grep (first column) and look if the numeric value