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