Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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_add, dbcsr_convert_dbcsr_to_csr, dbcsr_copy, dbcsr_create, &
17 : dbcsr_csr_create_from_dbcsr, dbcsr_csr_dbcsr_blkrow_dist, dbcsr_csr_destroy, &
18 : dbcsr_csr_type, dbcsr_csr_write, dbcsr_desymmetrize, dbcsr_finalize, dbcsr_get_block_p, &
19 : dbcsr_has_symmetry, dbcsr_p_type, dbcsr_put_block, dbcsr_release, dbcsr_set, dbcsr_type, &
20 : dbcsr_type_antisymmetric, dbcsr_type_complex_8, dbcsr_type_no_symmetry, dbcsr_type_real_8, &
21 : 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 16713 : 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 16713 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks
81 : TYPE(kpoint_type), POINTER :: kpoints
82 : TYPE(section_vals_type), POINTER :: dft_section
83 :
84 16713 : CALL timeset(routineN, handle)
85 :
86 : NULLIFY (dft_section)
87 :
88 16713 : logger => cp_get_default_logger()
89 16713 : output_unit = cp_logger_get_default_io_unit(logger)
90 :
91 16713 : 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 16713 : "PRINT%KS_CSR_WRITE"), cp_p_file)
94 :
95 16713 : 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 16713 : CALL timestop(handle)
110 :
111 16713 : 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 16713 : 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 16713 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
131 : TYPE(kpoint_type), POINTER :: kpoints
132 : TYPE(section_vals_type), POINTER :: dft_section
133 :
134 16713 : CALL timeset(routineN, handle)
135 :
136 : NULLIFY (dft_section)
137 :
138 16713 : logger => cp_get_default_logger()
139 16713 : output_unit = cp_logger_get_default_io_unit(logger)
140 :
141 16713 : 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 16713 : "PRINT%S_CSR_WRITE"), cp_p_file)
144 :
145 16713 : 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 16713 : CALL timestop(handle)
160 :
161 16713 : 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 : COMPLEX(KIND=dp), PARAMETER :: cone = CMPLX(1.0_dp, 0.0_dp, KIND=dp), &
324 : ione = CMPLX(0.0_dp, 1.0_dp, KIND=dp)
325 :
326 : CHARACTER(LEN=default_path_length) :: file_name, fileformat, subs_string
327 : INTEGER :: handle, igroup, ik, ikp, ispin, kplocal, &
328 : nkp_groups, output_unit, unit_nr
329 : INTEGER, DIMENSION(2) :: kp_range
330 20 : INTEGER, DIMENSION(:, :), POINTER :: kp_dist
331 : LOGICAL :: bin, uptr, use_real_wfn
332 : REAL(KIND=dp) :: thld
333 20 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
334 : TYPE(cp_logger_type), POINTER :: logger
335 : TYPE(dbcsr_csr_type) :: mat_csr
336 : TYPE(dbcsr_type) :: matrix_nosym
337 : TYPE(dbcsr_type), POINTER :: cmatrix, imatrix, imatrix_nosym, &
338 : rmatrix, rmatrix_nosym, tmatrix
339 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
340 20 : POINTER :: sab_nl
341 :
342 20 : CALL timeset(routineN, handle)
343 :
344 20 : logger => cp_get_default_logger()
345 20 : output_unit = cp_logger_get_default_io_unit(logger)
346 :
347 20 : subs_string = "PRINT%"//prefix//"_CSR_WRITE"
348 :
349 20 : CALL section_vals_val_get(dft_section, subs_string//"%THRESHOLD", r_val=thld)
350 20 : CALL section_vals_val_get(dft_section, subs_string//"%UPPER_TRIANGULAR", l_val=uptr)
351 20 : CALL section_vals_val_get(dft_section, subs_string//"%BINARY", l_val=bin)
352 :
353 20 : IF (bin) THEN
354 20 : fileformat = "UNFORMATTED"
355 : ELSE
356 0 : fileformat = "FORMATTED"
357 : END IF
358 :
359 20 : NULLIFY (sab_nl)
360 :
361 : ! Calculate the Hamiltonian at the k-points
362 : CALL get_kpoint_info(kpoints, xkp=xkp, use_real_wfn=use_real_wfn, kp_range=kp_range, &
363 20 : nkp_groups=nkp_groups, kp_dist=kp_dist, sab_nl=sab_nl)
364 :
365 20 : ALLOCATE (rmatrix)
366 : CALL dbcsr_create(rmatrix, template=mat(1, 1)%matrix, &
367 20 : matrix_type=dbcsr_type_symmetric)
368 20 : CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_nl)
369 :
370 20 : IF (.NOT. use_real_wfn) THEN
371 : ! Allocate temporary variables
372 16 : ALLOCATE (rmatrix_nosym, imatrix, imatrix_nosym, cmatrix, tmatrix)
373 : CALL dbcsr_create(rmatrix_nosym, template=mat(1, 1)%matrix, &
374 16 : matrix_type=dbcsr_type_no_symmetry)
375 : CALL dbcsr_create(imatrix, template=mat(1, 1)%matrix, &
376 16 : matrix_type=dbcsr_type_antisymmetric)
377 : CALL dbcsr_create(imatrix_nosym, template=mat(1, 1)%matrix, &
378 16 : matrix_type=dbcsr_type_no_symmetry)
379 : CALL dbcsr_create(cmatrix, template=mat(1, 1)%matrix, &
380 : matrix_type=dbcsr_type_no_symmetry, &
381 16 : data_type=dbcsr_type_complex_8)
382 : CALL dbcsr_create(tmatrix, template=mat(1, 1)%matrix, &
383 : matrix_type=dbcsr_type_no_symmetry, &
384 16 : data_type=dbcsr_type_complex_8)
385 16 : CALL cp_dbcsr_alloc_block_from_nbl(rmatrix_nosym, sab_nl)
386 16 : CALL cp_dbcsr_alloc_block_from_nbl(imatrix, sab_nl)
387 16 : CALL cp_dbcsr_alloc_block_from_nbl(imatrix_nosym, sab_nl)
388 16 : CALL cp_dbcsr_alloc_block_from_nbl(cmatrix, sab_nl)
389 16 : CALL cp_dbcsr_alloc_block_from_nbl(tmatrix, sab_nl)
390 : END IF
391 :
392 20 : kplocal = kp_range(2) - kp_range(1) + 1
393 56 : DO ikp = 1, kplocal
394 92 : DO ispin = 1, SIZE(mat, 1)
395 140 : DO igroup = 1, nkp_groups
396 : ! number of current kpoint
397 68 : ik = kp_dist(1, igroup) + ikp - 1
398 68 : CALL dbcsr_set(rmatrix, 0.0_dp)
399 68 : IF (use_real_wfn) THEN
400 : ! FT of the matrix
401 : CALL rskp_transform(rmatrix=rmatrix, rsmat=mat, ispin=ispin, &
402 4 : xkp=xkp(1:3, ik), cell_to_index=kpoints%cell_to_index, sab_nl=sab_nl)
403 : ! Convert to desymmetrized csr matrix
404 4 : CALL dbcsr_desymmetrize(rmatrix, matrix_nosym)
405 4 : CALL dbcsr_csr_create_from_dbcsr(matrix_nosym, mat_csr, dbcsr_csr_dbcsr_blkrow_dist)
406 4 : CALL dbcsr_convert_dbcsr_to_csr(matrix_nosym, mat_csr)
407 4 : CALL dbcsr_release(matrix_nosym)
408 : ELSE
409 : ! FT of the matrix
410 64 : CALL dbcsr_set(imatrix, 0.0_dp)
411 : CALL rskp_transform(rmatrix=rmatrix, cmatrix=imatrix, rsmat=mat, ispin=ispin, &
412 64 : xkp=xkp(1:3, ik), cell_to_index=kpoints%cell_to_index, sab_nl=sab_nl)
413 :
414 : ! Desymmetrize and sum the real and imaginary part into
415 : ! cmatrix
416 64 : CALL dbcsr_desymmetrize(rmatrix, rmatrix_nosym)
417 64 : CALL dbcsr_desymmetrize(imatrix, imatrix_nosym)
418 64 : CALL dbcsr_copy(cmatrix, rmatrix_nosym)
419 64 : CALL dbcsr_copy(tmatrix, imatrix_nosym)
420 64 : CALL dbcsr_add(cmatrix, tmatrix, cone, ione)
421 : ! Convert to csr
422 64 : CALL dbcsr_csr_create_from_dbcsr(cmatrix, mat_csr, dbcsr_csr_dbcsr_blkrow_dist)
423 64 : CALL dbcsr_convert_dbcsr_to_csr(cmatrix, mat_csr)
424 : END IF
425 : ! Write to file
426 68 : WRITE (file_name, '(2(A,I0))') prefix//"_SPIN_", ispin, "_K_", ik
427 : unit_nr = cp_print_key_unit_nr(logger, dft_section, subs_string, &
428 : extension=".csr", middle_name=TRIM(file_name), &
429 68 : file_status="REPLACE", file_form=fileformat)
430 68 : CALL dbcsr_csr_write(mat_csr, unit_nr, upper_triangle=uptr, threshold=thld, binary=bin)
431 :
432 68 : CALL cp_print_key_finished_output(unit_nr, logger, dft_section, subs_string)
433 :
434 104 : CALL dbcsr_csr_destroy(mat_csr)
435 : END DO
436 : END DO
437 : END DO
438 20 : CALL dbcsr_release(rmatrix)
439 20 : DEALLOCATE (rmatrix)
440 20 : IF (.NOT. use_real_wfn) THEN
441 16 : CALL dbcsr_release(rmatrix_nosym)
442 16 : CALL dbcsr_release(imatrix)
443 16 : CALL dbcsr_release(imatrix_nosym)
444 16 : CALL dbcsr_release(cmatrix)
445 16 : CALL dbcsr_release(tmatrix)
446 16 : DEALLOCATE (rmatrix_nosym, imatrix, imatrix_nosym, cmatrix, tmatrix)
447 : END IF
448 20 : CALL timestop(handle)
449 :
450 20 : END SUBROUTINE write_matrix_kp_csr
451 :
452 : ! **************************************************************************************************
453 : !> \brief Desymmetrizes the KS or S matrices which are stored in symmetric !matrices
454 : !> \param mat Hamiltonian or overlap matrices
455 : !> \param mat_nosym Desymmetrized Hamiltonian or overlap matrices
456 : !> \param cell_to_index Mapping of cell indices to linear RS indices
457 : !> \param index_to_cell Mapping of linear RS indices to cell indices
458 : !> \param kpoints Kpoint environment
459 : !> \author Fabian Ducry
460 : ! **************************************************************************************************
461 2 : SUBROUTINE desymmetrize_rs_matrix(mat, mat_nosym, cell_to_index, index_to_cell, kpoints)
462 : TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(IN), &
463 : POINTER :: mat
464 : TYPE(dbcsr_p_type), DIMENSION(:, :), &
465 : INTENT(INOUT), POINTER :: mat_nosym
466 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :), &
467 : INTENT(OUT) :: cell_to_index
468 : INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(OUT) :: index_to_cell
469 : TYPE(kpoint_type), INTENT(IN), POINTER :: kpoints
470 :
471 : CHARACTER(len=*), PARAMETER :: routineN = 'desymmetrize_rs_matrix'
472 :
473 : INTEGER :: handle, iatom, ic, icn, icol, irow, &
474 : ispin, jatom, ncell, nomirror, nspin, &
475 : nx, ny, nz
476 : INTEGER, DIMENSION(3) :: cell
477 2 : INTEGER, DIMENSION(:, :), POINTER :: i2c
478 2 : INTEGER, DIMENSION(:, :, :), POINTER :: c2i
479 : LOGICAL :: found, lwtr
480 2 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
481 : TYPE(neighbor_list_iterator_p_type), &
482 2 : DIMENSION(:), POINTER :: nl_iterator
483 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
484 2 : POINTER :: sab_nl
485 :
486 2 : CALL timeset(routineN, handle)
487 :
488 2 : i2c => kpoints%index_to_cell
489 2 : c2i => kpoints%cell_to_index
490 :
491 2 : ncell = SIZE(i2c, 2)
492 2 : nspin = SIZE(mat, 1)
493 :
494 2 : nx = MAX(ABS(LBOUND(c2i, 1)), ABS(UBOUND(c2i, 1)))
495 2 : ny = MAX(ABS(LBOUND(c2i, 2)), ABS(UBOUND(c2i, 3)))
496 2 : nz = MAX(ABS(LBOUND(c2i, 3)), ABS(UBOUND(c2i, 3)))
497 10 : ALLOCATE (cell_to_index(-nx:nx, -ny:ny, -nz:nz))
498 : cell_to_index(LBOUND(c2i, 1):UBOUND(c2i, 1), &
499 : LBOUND(c2i, 2):UBOUND(c2i, 2), &
500 84 : LBOUND(c2i, 3):UBOUND(c2i, 3)) = c2i
501 :
502 : ! identify cells with no mirror img
503 : nomirror = 0
504 52 : DO ic = 1, ncell
505 200 : cell = i2c(:, ic)
506 50 : IF (cell_to_index(-cell(1), -cell(2), -cell(3)) == 0) &
507 6 : nomirror = nomirror + 1
508 : END DO
509 :
510 : ! create the mirror imgs
511 6 : ALLOCATE (index_to_cell(3, ncell + nomirror))
512 202 : index_to_cell(:, 1:ncell) = i2c
513 :
514 : nomirror = 0 ! count the imgs without mirror
515 52 : DO ic = 1, ncell
516 200 : cell = index_to_cell(:, ic)
517 52 : IF (cell_to_index(-cell(1), -cell(2), -cell(3)) == 0) THEN
518 4 : nomirror = nomirror + 1
519 16 : index_to_cell(:, ncell + nomirror) = -cell
520 4 : cell_to_index(-cell(1), -cell(2), -cell(3)) = ncell + nomirror
521 : END IF
522 : END DO
523 2 : ncell = ncell + nomirror
524 :
525 2 : CALL get_kpoint_info(kpoints, sab_nl=sab_nl)
526 : ! allocate the nonsymmetric matrices
527 2 : NULLIFY (mat_nosym)
528 2 : CALL dbcsr_allocate_matrix_set(mat_nosym, nspin, ncell)
529 4 : DO ispin = 1, nspin
530 58 : DO ic = 1, ncell
531 54 : ALLOCATE (mat_nosym(ispin, ic)%matrix)
532 : CALL dbcsr_create(matrix=mat_nosym(ispin, ic)%matrix, &
533 : template=mat(1, 1)%matrix, &
534 : matrix_type=dbcsr_type_no_symmetry, &
535 54 : data_type=dbcsr_type_real_8)
536 : CALL cp_dbcsr_alloc_block_from_nbl(mat_nosym(ispin, ic)%matrix, &
537 54 : sab_nl, desymmetrize=.TRUE.)
538 56 : CALL dbcsr_set(mat_nosym(ispin, ic)%matrix, 0.0_dp)
539 : END DO
540 : END DO
541 :
542 4 : DO ispin = 1, nspin
543 : ! desymmetrize the matrix for real space printing
544 2 : CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
545 506 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
546 504 : CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom, cell=cell)
547 :
548 504 : ic = cell_to_index(cell(1), cell(2), cell(3))
549 504 : icn = cell_to_index(-cell(1), -cell(2), -cell(3))
550 504 : CPASSERT(icn > 0)
551 :
552 504 : irow = iatom
553 504 : icol = jatom
554 504 : lwtr = .FALSE.
555 : ! always copy from the top
556 504 : IF (iatom > jatom) THEN
557 216 : irow = jatom
558 216 : icol = iatom
559 216 : lwtr = .TRUE.
560 : END IF
561 :
562 : CALL dbcsr_get_block_p(matrix=mat(ispin, ic)%matrix, &
563 504 : row=irow, col=icol, block=block, found=found)
564 504 : CPASSERT(found)
565 :
566 : ! copy to M(R) at (iatom,jatom)
567 : ! copy to M(-R) at (jatom,iatom)
568 506 : IF (lwtr) THEN
569 : CALL dbcsr_put_block(matrix=mat_nosym(ispin, ic)%matrix, &
570 4536 : row=iatom, col=jatom, block=TRANSPOSE(block))
571 : CALL dbcsr_put_block(matrix=mat_nosym(ispin, icn)%matrix, &
572 216 : row=jatom, col=iatom, block=block)
573 : ELSE
574 : CALL dbcsr_put_block(matrix=mat_nosym(ispin, ic)%matrix, &
575 288 : row=iatom, col=jatom, block=block)
576 : CALL dbcsr_put_block(matrix=mat_nosym(ispin, icn)%matrix, &
577 6048 : row=jatom, col=iatom, block=TRANSPOSE(block))
578 : END IF
579 : END DO
580 4 : CALL neighbor_list_iterator_release(nl_iterator)
581 : END DO
582 :
583 4 : DO ispin = 1, nspin
584 58 : DO ic = 1, ncell
585 56 : CALL dbcsr_finalize(mat_nosym(ispin, ic)%matrix)
586 : END DO
587 : END DO
588 :
589 2 : CALL timestop(handle)
590 :
591 2 : END SUBROUTINE desymmetrize_rs_matrix
592 :
593 : END MODULE qs_scf_csr_write
|
