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