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 Restart file for k point calculations
10 : !> \par History
11 : !> \author JGH (30.09.2015)
12 : ! **************************************************************************************************
13 : MODULE kpoint_io
14 :
15 : USE atomic_kind_types, ONLY: get_atomic_kind
16 : USE basis_set_types, ONLY: get_gto_basis_set,&
17 : gto_basis_set_type
18 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
19 : copy_fm_to_dbcsr
20 : USE cp_files, ONLY: close_file,&
21 : open_file
22 : USE cp_fm_types, ONLY: cp_fm_read_unformatted,&
23 : cp_fm_type,&
24 : cp_fm_write_unformatted
25 : USE cp_log_handling, ONLY: cp_get_default_logger,&
26 : cp_logger_type,&
27 : cp_to_string
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 dbcsr_api, ONLY: dbcsr_get_info,&
33 : dbcsr_p_type
34 : USE input_section_types, ONLY: section_vals_type
35 : USE kinds, ONLY: default_path_length
36 : USE kpoint_types, ONLY: get_kpoint_info,&
37 : kpoint_type
38 : USE message_passing, ONLY: mp_para_env_type
39 : USE orbital_pointers, ONLY: nso
40 : USE particle_types, ONLY: particle_type
41 : USE qs_dftb_types, ONLY: qs_dftb_atom_type
42 : USE qs_dftb_utils, ONLY: get_dftb_atom_param
43 : USE qs_kind_types, ONLY: get_qs_kind,&
44 : qs_kind_type
45 : USE qs_mo_io, ONLY: wfn_restart_file_name
46 : USE qs_scf_types, ONLY: qs_scf_env_type
47 : #include "./base/base_uses.f90"
48 :
49 : IMPLICIT NONE
50 :
51 : PRIVATE
52 :
53 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'kpoint_io'
54 :
55 : INTEGER, PARAMETER :: version = 1
56 :
57 : PUBLIC :: read_kpoints_restart, &
58 : write_kpoints_restart
59 :
60 : ! **************************************************************************************************
61 :
62 : CONTAINS
63 :
64 : ! **************************************************************************************************
65 : !> \brief ...
66 : !> \param denmat ...
67 : !> \param kpoints ...
68 : !> \param scf_env ...
69 : !> \param dft_section ...
70 : !> \param particle_set ...
71 : !> \param qs_kind_set ...
72 : ! **************************************************************************************************
73 5320 : SUBROUTINE write_kpoints_restart(denmat, kpoints, scf_env, dft_section, &
74 : particle_set, qs_kind_set)
75 :
76 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: denmat
77 : TYPE(kpoint_type), POINTER :: kpoints
78 : TYPE(qs_scf_env_type), POINTER :: scf_env
79 : TYPE(section_vals_type), POINTER :: dft_section
80 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
81 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
82 :
83 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_kpoints_restart'
84 : CHARACTER(LEN=30), DIMENSION(2), PARAMETER :: &
85 : keys = (/"SCF%PRINT%RESTART_HISTORY", "SCF%PRINT%RESTART "/)
86 :
87 : INTEGER :: handle, ic, ikey, ires, ispin, nimages, &
88 : nspin
89 : INTEGER, DIMENSION(3) :: cell
90 5320 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
91 : TYPE(cp_fm_type), POINTER :: fmat
92 : TYPE(cp_logger_type), POINTER :: logger
93 :
94 5320 : CALL timeset(routineN, handle)
95 5320 : logger => cp_get_default_logger()
96 :
97 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
98 5320 : dft_section, keys(1)), cp_p_file) .OR. &
99 : BTEST(cp_print_key_should_output(logger%iter_info, &
100 : dft_section, keys(2)), cp_p_file)) THEN
101 :
102 1434 : DO ikey = 1, SIZE(keys)
103 :
104 956 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
105 5320 : dft_section, keys(ikey)), cp_p_file)) THEN
106 :
107 : ires = cp_print_key_unit_nr(logger, dft_section, keys(ikey), &
108 : extension=".kp", file_status="REPLACE", file_action="WRITE", &
109 478 : do_backup=.TRUE., file_form="UNFORMATTED")
110 :
111 478 : CALL write_kpoints_file_header(qs_kind_set, particle_set, ires)
112 :
113 478 : nspin = SIZE(denmat, 1)
114 478 : nimages = SIZE(denmat, 2)
115 478 : NULLIFY (cell_to_index)
116 478 : IF (nimages > 1) THEN
117 478 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
118 : END IF
119 478 : CPASSERT(ASSOCIATED(scf_env%scf_work1))
120 478 : fmat => scf_env%scf_work1(1)
121 :
122 1070 : DO ispin = 1, nspin
123 592 : IF (ires > 0) WRITE (ires) ispin, nspin, nimages
124 50618 : DO ic = 1, nimages
125 49548 : IF (nimages > 1) THEN
126 49548 : cell = get_cell(ic, cell_to_index)
127 : ELSE
128 0 : cell = 0
129 : END IF
130 49548 : IF (ires > 0) WRITE (ires) ic, cell
131 49548 : CALL copy_dbcsr_to_fm(denmat(ispin, ic)%matrix, fmat)
132 50140 : CALL cp_fm_write_unformatted(fmat, ires)
133 : END DO
134 : END DO
135 :
136 478 : CALL cp_print_key_finished_output(ires, logger, dft_section, TRIM(keys(ikey)))
137 :
138 : END IF
139 :
140 : END DO
141 :
142 : END IF
143 :
144 5320 : CALL timestop(handle)
145 :
146 5320 : END SUBROUTINE write_kpoints_restart
147 :
148 : ! **************************************************************************************************
149 : !> \brief ...
150 : !> \param ic ...
151 : !> \param cell_to_index ...
152 : !> \return ...
153 : ! **************************************************************************************************
154 49548 : FUNCTION get_cell(ic, cell_to_index) RESULT(cell)
155 : INTEGER, INTENT(IN) :: ic
156 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
157 : INTEGER, DIMENSION(3) :: cell
158 :
159 : INTEGER :: i1, i2, i3
160 :
161 198192 : cell = 0
162 257728 : ALLCELL: DO i3 = LBOUND(cell_to_index, 3), UBOUND(cell_to_index, 3)
163 1508146 : DO i2 = LBOUND(cell_to_index, 2), UBOUND(cell_to_index, 2)
164 11903314 : DO i1 = LBOUND(cell_to_index, 1), UBOUND(cell_to_index, 1)
165 9630634 : IF (cell_to_index(i1, i2, i3) == ic) THEN
166 49548 : cell(1) = i1
167 49548 : cell(2) = i2
168 49548 : cell(3) = i3
169 49548 : EXIT ALLCELL
170 : END IF
171 : END DO
172 : END DO
173 : END DO ALLCELL
174 :
175 49548 : END FUNCTION get_cell
176 :
177 : ! **************************************************************************************************
178 : !> \brief ...
179 : !> \param qs_kind_set ...
180 : !> \param particle_set ...
181 : !> \param ires ...
182 : ! **************************************************************************************************
183 478 : SUBROUTINE write_kpoints_file_header(qs_kind_set, particle_set, ires)
184 :
185 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
186 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
187 : INTEGER :: ires
188 :
189 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_kpoints_file_header'
190 :
191 : INTEGER :: handle, iatom, ikind, iset, ishell, &
192 : lmax, lshell, nao, natom, nset, &
193 : nset_max, nsgf, nshell_max
194 478 : INTEGER, DIMENSION(:), POINTER :: nset_info, nshell
195 478 : INTEGER, DIMENSION(:, :), POINTER :: l, nshell_info
196 478 : INTEGER, DIMENSION(:, :, :), POINTER :: nso_info
197 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
198 : TYPE(qs_dftb_atom_type), POINTER :: dftb_parameter
199 :
200 478 : CALL timeset(routineN, handle)
201 :
202 478 : IF (ires > 0) THEN
203 : ! create some info about the basis set first
204 239 : natom = SIZE(particle_set, 1)
205 239 : nset_max = 0
206 239 : nshell_max = 0
207 :
208 1736 : DO iatom = 1, natom
209 1497 : NULLIFY (orb_basis_set, dftb_parameter)
210 1497 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
211 : CALL get_qs_kind(qs_kind_set(ikind), &
212 1497 : basis_set=orb_basis_set, dftb_parameter=dftb_parameter)
213 3233 : IF (ASSOCIATED(orb_basis_set)) THEN
214 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
215 : nset=nset, &
216 : nshell=nshell, &
217 1041 : l=l)
218 1041 : nset_max = MAX(nset_max, nset)
219 3302 : DO iset = 1, nset
220 3302 : nshell_max = MAX(nshell_max, nshell(iset))
221 : END DO
222 456 : ELSEIF (ASSOCIATED(dftb_parameter)) THEN
223 456 : CALL get_dftb_atom_param(dftb_parameter, lmax=lmax)
224 456 : nset_max = MAX(nset_max, 1)
225 456 : nshell_max = MAX(nshell_max, lmax + 1)
226 : ELSE
227 0 : CPABORT("Unknown basis type.")
228 : END IF
229 : END DO
230 :
231 1195 : ALLOCATE (nso_info(nshell_max, nset_max, natom))
232 8276 : nso_info(:, :, :) = 0
233 :
234 956 : ALLOCATE (nshell_info(nset_max, natom))
235 4460 : nshell_info(:, :) = 0
236 :
237 717 : ALLOCATE (nset_info(natom))
238 1736 : nset_info(:) = 0
239 :
240 239 : nao = 0
241 1736 : DO iatom = 1, natom
242 1497 : NULLIFY (orb_basis_set, dftb_parameter)
243 1497 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
244 : CALL get_qs_kind(qs_kind_set(ikind), &
245 1497 : basis_set=orb_basis_set, dftb_parameter=dftb_parameter)
246 3233 : IF (ASSOCIATED(orb_basis_set)) THEN
247 1041 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, nsgf=nsgf, nset=nset, nshell=nshell, l=l)
248 1041 : nao = nao + nsgf
249 1041 : nset_info(iatom) = nset
250 3302 : DO iset = 1, nset
251 2261 : nshell_info(iset, iatom) = nshell(iset)
252 5681 : DO ishell = 1, nshell(iset)
253 2379 : lshell = l(ishell, iset)
254 4640 : nso_info(ishell, iset, iatom) = nso(lshell)
255 : END DO
256 : END DO
257 456 : ELSEIF (ASSOCIATED(dftb_parameter)) THEN
258 456 : CALL get_dftb_atom_param(dftb_parameter, lmax=lmax)
259 456 : nset_info(iatom) = 1
260 456 : nshell_info(1, iatom) = lmax + 1
261 1768 : DO ishell = 1, lmax + 1
262 1312 : lshell = ishell - 1
263 1768 : nso_info(ishell, 1, iatom) = nso(lshell)
264 : END DO
265 456 : nao = nao + (lmax + 1)**2
266 : ELSE
267 0 : CPABORT("Unknown basis set type.")
268 : END IF
269 : END DO
270 :
271 239 : WRITE (ires) version
272 239 : WRITE (ires) natom, nao, nset_max, nshell_max
273 1736 : WRITE (ires) nset_info
274 4460 : WRITE (ires) nshell_info
275 8276 : WRITE (ires) nso_info
276 :
277 239 : DEALLOCATE (nset_info, nshell_info, nso_info)
278 : END IF
279 :
280 478 : CALL timestop(handle)
281 :
282 478 : END SUBROUTINE write_kpoints_file_header
283 :
284 : ! **************************************************************************************************
285 : !> \brief ...
286 : !> \param denmat ...
287 : !> \param kpoints ...
288 : !> \param fmwork ...
289 : !> \param natom ...
290 : !> \param para_env ...
291 : !> \param id_nr ...
292 : !> \param dft_section ...
293 : !> \param natom_mismatch ...
294 : ! **************************************************************************************************
295 6 : SUBROUTINE read_kpoints_restart(denmat, kpoints, fmwork, natom, &
296 : para_env, id_nr, dft_section, natom_mismatch)
297 :
298 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: denmat
299 : TYPE(kpoint_type), POINTER :: kpoints
300 : TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: fmwork
301 : INTEGER, INTENT(IN) :: natom
302 : TYPE(mp_para_env_type), POINTER :: para_env
303 : INTEGER, INTENT(IN) :: id_nr
304 : TYPE(section_vals_type), POINTER :: dft_section
305 : LOGICAL, INTENT(OUT) :: natom_mismatch
306 :
307 : CHARACTER(LEN=*), PARAMETER :: routineN = 'read_kpoints_restart'
308 :
309 : CHARACTER(LEN=default_path_length) :: file_name
310 : INTEGER :: handle, restart_unit
311 : LOGICAL :: exist
312 : TYPE(cp_logger_type), POINTER :: logger
313 :
314 6 : CALL timeset(routineN, handle)
315 6 : logger => cp_get_default_logger()
316 :
317 6 : restart_unit = -1
318 :
319 6 : IF (para_env%is_source()) THEN
320 :
321 3 : CALL wfn_restart_file_name(file_name, exist, dft_section, logger, kp=.TRUE.)
322 3 : IF (id_nr /= 0) THEN
323 : ! Is it one of the backup files?
324 0 : file_name = TRIM(file_name)//".bak-"//ADJUSTL(cp_to_string(id_nr))
325 : END IF
326 :
327 : CALL open_file(file_name=file_name, &
328 : file_action="READ", &
329 : file_form="UNFORMATTED", &
330 : file_status="OLD", &
331 3 : unit_number=restart_unit)
332 :
333 : END IF
334 :
335 : CALL read_kpoints_restart_low(denmat, kpoints, fmwork(1), para_env, &
336 6 : natom, restart_unit, natom_mismatch)
337 :
338 : ! only the io_node returns natom_mismatch, must broadcast it
339 6 : CALL para_env%bcast(natom_mismatch)
340 :
341 : ! Close restart file
342 6 : IF (para_env%is_source()) CALL close_file(unit_number=restart_unit)
343 :
344 6 : CALL timestop(handle)
345 :
346 6 : END SUBROUTINE read_kpoints_restart
347 :
348 : ! **************************************************************************************************
349 : !> \brief Reading the mos from apreviously defined restart file
350 : !> \param denmat ...
351 : !> \param kpoints ...
352 : !> \param fmat ...
353 : !> \param para_env ...
354 : !> \param natom ...
355 : !> \param rst_unit ...
356 : !> \param natom_mismatch ...
357 : !> \par History
358 : !> 12.2007 created [MI]
359 : !> \author MI
360 : ! **************************************************************************************************
361 12 : SUBROUTINE read_kpoints_restart_low(denmat, kpoints, fmat, para_env, &
362 : natom, rst_unit, natom_mismatch)
363 :
364 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: denmat
365 : TYPE(kpoint_type), POINTER :: kpoints
366 : TYPE(cp_fm_type), INTENT(INOUT) :: fmat
367 : TYPE(mp_para_env_type), POINTER :: para_env
368 : INTEGER, INTENT(IN) :: natom, rst_unit
369 : LOGICAL, INTENT(OUT) :: natom_mismatch
370 :
371 : INTEGER :: ic, image, ispin, nao, nao_read, &
372 : natom_read, ni, nimages, nset_max, &
373 : nshell_max, nspin, nspin_read, &
374 : version_read
375 : INTEGER, DIMENSION(3) :: cell
376 6 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
377 : LOGICAL :: natom_match
378 :
379 0 : CPASSERT(ASSOCIATED(denmat))
380 6 : CALL dbcsr_get_info(denmat(1, 1)%matrix, nfullrows_total=nao)
381 :
382 6 : nspin = SIZE(denmat, 1)
383 6 : nimages = SIZE(denmat, 2)
384 6 : NULLIFY (cell_to_index)
385 6 : IF (nimages > 1) THEN
386 6 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
387 : END IF
388 :
389 6 : IF (para_env%is_source()) THEN
390 3 : READ (rst_unit) version_read
391 3 : IF (version_read /= version) THEN
392 : CALL cp_abort(__LOCATION__, &
393 0 : " READ RESTART : Version of restart file not supported")
394 : END IF
395 3 : READ (rst_unit) natom_read, nao_read, nset_max, nshell_max
396 3 : natom_match = (natom_read == natom)
397 3 : IF (natom_match) THEN
398 3 : IF (nao_read /= nao) THEN
399 : CALL cp_abort(__LOCATION__, &
400 0 : " READ RESTART : Different number of basis functions")
401 : END IF
402 3 : READ (rst_unit)
403 3 : READ (rst_unit)
404 3 : READ (rst_unit)
405 : END IF
406 : END IF
407 :
408 : ! make natom_match and natom_mismatch uniform across all nodes
409 6 : CALL para_env%bcast(natom_match)
410 6 : natom_mismatch = .NOT. natom_match
411 : ! handle natom_match false
412 6 : IF (natom_mismatch) THEN
413 : CALL cp_warn(__LOCATION__, &
414 0 : " READ RESTART : WARNING : DIFFERENT natom, returning ")
415 : ELSE
416 :
417 : DO
418 10 : ispin = 0
419 10 : IF (para_env%is_source()) THEN
420 5 : READ (rst_unit) ispin, nspin_read, ni
421 : END IF
422 10 : CALL para_env%bcast(ispin)
423 10 : CALL para_env%bcast(nspin_read)
424 10 : CALL para_env%bcast(ni)
425 10 : IF (nspin_read /= nspin) THEN
426 : CALL cp_abort(__LOCATION__, &
427 0 : " READ RESTART : Different spin polarisation not supported")
428 : END IF
429 816 : DO ic = 1, ni
430 806 : IF (para_env%is_source()) THEN
431 403 : READ (rst_unit) image, cell
432 : END IF
433 806 : CALL para_env%bcast(image)
434 806 : CALL para_env%bcast(cell)
435 : !
436 806 : CALL cp_fm_read_unformatted(fmat, rst_unit)
437 : !
438 806 : IF (nimages > 1) THEN
439 806 : image = get_index(cell, cell_to_index)
440 0 : ELSE IF (SUM(ABS(cell(:))) == 0) THEN
441 0 : image = 1
442 : ELSE
443 0 : image = 0
444 : END IF
445 816 : IF (image >= 1 .AND. image <= nimages) THEN
446 806 : CALL copy_fm_to_dbcsr(fmat, denmat(ispin, image)%matrix)
447 : END IF
448 : END DO
449 10 : IF (ispin == nspin_read) EXIT
450 : END DO
451 :
452 : END IF
453 :
454 6 : END SUBROUTINE read_kpoints_restart_low
455 :
456 : ! **************************************************************************************************
457 : !> \brief ...
458 : !> \param cell ...
459 : !> \param cell_to_index ...
460 : !> \return ...
461 : ! **************************************************************************************************
462 806 : FUNCTION get_index(cell, cell_to_index) RESULT(index)
463 : INTEGER, DIMENSION(3), INTENT(IN) :: cell
464 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
465 : INTEGER :: index
466 :
467 : IF (CELL(1) >= LBOUND(cell_to_index, 1) .AND. CELL(1) <= UBOUND(cell_to_index, 1) .AND. &
468 : CELL(2) >= LBOUND(cell_to_index, 2) .AND. CELL(2) <= UBOUND(cell_to_index, 2) .AND. &
469 5642 : CELL(3) >= LBOUND(cell_to_index, 3) .AND. CELL(3) <= UBOUND(cell_to_index, 3)) THEN
470 :
471 806 : index = cell_to_index(cell(1), cell(2), cell(3))
472 :
473 : ELSE
474 :
475 : index = 0
476 :
477 : END IF
478 :
479 806 : END FUNCTION get_index
480 :
481 : END MODULE kpoint_io
|
