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_real_8, 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 17161 : 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 17161 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks
81 : TYPE(kpoint_type), POINTER :: kpoints
82 : TYPE(section_vals_type), POINTER :: dft_section
83 :
84 17161 : CALL timeset(routineN, handle)
85 :
86 : NULLIFY (dft_section)
87 :
88 17161 : logger => cp_get_default_logger()
89 17161 : output_unit = cp_logger_get_default_io_unit(logger)
90 :
91 17161 : 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 17161 : "PRINT%KS_CSR_WRITE"), cp_p_file)
94 :
95 17161 : 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 17161 : CALL timestop(handle)
110 :
111 17161 : 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 17161 : 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 17161 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
131 : TYPE(kpoint_type), POINTER :: kpoints
132 : TYPE(section_vals_type), POINTER :: dft_section
133 :
134 17161 : CALL timeset(routineN, handle)
135 :
136 : NULLIFY (dft_section)
137 :
138 17161 : logger => cp_get_default_logger()
139 17161 : output_unit = cp_logger_get_default_io_unit(logger)
140 :
141 17161 : 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 17161 : "PRINT%S_CSR_WRITE"), cp_p_file)
144 :
145 17161 : 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 17161 : CALL timestop(handle)
160 :
161 17161 : 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 : matrix_type=dbcsr_type_no_symmetry, &
522 54 : data_type=dbcsr_type_real_8)
523 : CALL cp_dbcsr_alloc_block_from_nbl(mat_nosym(ispin, ic)%matrix, &
524 54 : sab_nl, desymmetrize=.TRUE.)
525 56 : CALL dbcsr_set(mat_nosym(ispin, ic)%matrix, 0.0_dp)
526 : END DO
527 : END DO
528 :
529 4 : DO ispin = 1, nspin
530 : ! desymmetrize the matrix for real space printing
531 2 : CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
532 506 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
533 504 : CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom, cell=cell)
534 :
535 504 : ic = cell_to_index(cell(1), cell(2), cell(3))
536 504 : icn = cell_to_index(-cell(1), -cell(2), -cell(3))
537 504 : CPASSERT(icn > 0)
538 :
539 504 : irow = iatom
540 504 : icol = jatom
541 504 : lwtr = .FALSE.
542 : ! always copy from the top
543 504 : IF (iatom > jatom) THEN
544 216 : irow = jatom
545 216 : icol = iatom
546 216 : lwtr = .TRUE.
547 : END IF
548 :
549 : CALL dbcsr_get_block_p(matrix=mat(ispin, ic)%matrix, &
550 504 : row=irow, col=icol, block=block, found=found)
551 504 : CPASSERT(found)
552 :
553 : ! copy to M(R) at (iatom,jatom)
554 : ! copy to M(-R) at (jatom,iatom)
555 506 : IF (lwtr) THEN
556 : CALL dbcsr_put_block(matrix=mat_nosym(ispin, ic)%matrix, &
557 4536 : row=iatom, col=jatom, block=TRANSPOSE(block))
558 : CALL dbcsr_put_block(matrix=mat_nosym(ispin, icn)%matrix, &
559 216 : row=jatom, col=iatom, block=block)
560 : ELSE
561 : CALL dbcsr_put_block(matrix=mat_nosym(ispin, ic)%matrix, &
562 288 : row=iatom, col=jatom, block=block)
563 : CALL dbcsr_put_block(matrix=mat_nosym(ispin, icn)%matrix, &
564 6048 : row=jatom, col=iatom, block=TRANSPOSE(block))
565 : END IF
566 : END DO
567 4 : CALL neighbor_list_iterator_release(nl_iterator)
568 : END DO
569 :
570 4 : DO ispin = 1, nspin
571 58 : DO ic = 1, ncell
572 56 : CALL dbcsr_finalize(mat_nosym(ispin, ic)%matrix)
573 : END DO
574 : END DO
575 :
576 2 : CALL timestop(handle)
577 :
578 2 : END SUBROUTINE desymmetrize_rs_matrix
579 :
580 : END MODULE qs_scf_csr_write
|
