1!! Copyright (C) 2002-2016 M. Marques, A. Castro, A. Rubio, G. Bertsch, X. Andrade 2!! 3!! This program is free software; you can redistribute it and/or modify 4!! it under the terms of the GNU General Public License as published by 5!! the Free Software Foundation; either version 2, or (at your option) 6!! any later version. 7!! 8!! This program is distributed in the hope that it will be useful, 9!! but WITHOUT ANY WARRANTY; without even the implied warranty of 10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11!! GNU General Public License for more details. 12!! 13!! You should have received a copy of the GNU General Public License 14!! along with this program; if not, write to the Free Software 15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 16!! 02110-1301, USA. 17!! 18 19#include "global.h" 20 21module states_elec_restart_oct_m 22 use global_oct_m 23 use grid_oct_m 24 use io_function_oct_m 25 use kpoints_oct_m 26 use lalg_basic_oct_m 27 use linear_response_oct_m 28 use loct_oct_m 29 use mesh_oct_m 30 use mesh_function_oct_m 31 use messages_oct_m 32 use mpi_oct_m 33 use multicomm_oct_m 34 use multigrid_oct_m 35 use namespace_oct_m 36 use parser_oct_m 37 use profiling_oct_m 38 use restart_oct_m 39 use simul_box_oct_m 40 use smear_oct_m 41 use states_abst_oct_m 42 use states_elec_oct_m 43 use states_elec_dim_oct_m 44 use string_oct_m 45 use types_oct_m 46 47 implicit none 48 49 type(profile_t), save :: prof_read, prof_write 50 51 private 52 public :: & 53 states_elec_look_and_load, & 54 states_elec_dump, & 55 states_elec_load, & 56 states_elec_dump_rho, & 57 states_elec_load_rho, & 58 states_elec_dump_spin, & 59 states_elec_load_spin, & 60 states_elec_dump_frozen, & 61 states_elec_load_frozen, & 62 states_elec_read_user_def_orbitals 63 64contains 65 66 ! --------------------------------------------------------- 67 subroutine states_elec_look_and_load(restart, namespace, st, gr, is_complex) 68 type(restart_t), intent(in) :: restart 69 type(namespace_t), intent(in) :: namespace 70 type(states_elec_t), target, intent(inout) :: st 71 type(grid_t), intent(in) :: gr 72 logical, optional, intent(in) :: is_complex 73 74 integer :: kpoints, dim, nst, ierr 75 FLOAT, pointer :: new_occ(:,:) 76 77 PUSH_SUB(states_elec_look_and_load) 78 79 !check how many wfs we have 80 call states_elec_look(restart, kpoints, dim, nst, ierr) 81 if(ierr /= 0) then 82 message(1) = "Unable to read states information." 83 call messages_fatal(1, namespace=namespace) 84 end if 85 86 if(st%parallel_in_states) then 87 message(1) = "Internal error: cannot use states_elec_look_and_load when parallel in states." 88 call messages_fatal(1, namespace=namespace) 89 end if 90 91 ! Resize st%occ, retaining information there 92 SAFE_ALLOCATE(new_occ(1:nst, 1:st%d%nik)) 93 new_occ(:,:) = M_ZERO 94 new_occ(1:min(nst, st%nst),:) = st%occ(1:min(nst, st%nst),:) 95 SAFE_DEALLOCATE_P(st%occ) 96 st%occ => new_occ 97 98 ! FIXME: This wrong, one cannot just change the number of states 99 ! without updating the internal structures, in the case of parallel in states. 100 ! I think it is right now without state parallelism. 101 st%nst = nst 102 st%st_start = 1 103 st%st_end = nst 104 st%lnst = nst 105 106 SAFE_DEALLOCATE_P(st%node) 107 SAFE_ALLOCATE(st%node(1:st%nst)) 108 st%node(:) = 0 109 110 SAFE_DEALLOCATE_P(st%eigenval) 111 SAFE_ALLOCATE(st%eigenval(1:st%nst, 1:st%d%nik)) 112 st%eigenval = huge(st%eigenval) 113 114 if (present(is_complex)) then 115 if ( is_complex ) then 116 call states_elec_allocate_wfns(st, gr%mesh, TYPE_CMPLX) 117 else 118 call states_elec_allocate_wfns(st, gr%mesh, TYPE_FLOAT) 119 end if 120 else 121 ! allow states_elec_allocate_wfns to decide for itself whether complex or real needed 122 call states_elec_allocate_wfns(st, gr%mesh) 123 end if 124 125 if(st%d%ispin == SPINORS) then 126 SAFE_ALLOCATE(st%spin(1:3, 1:st%nst, 1:st%d%nik)) 127 st%spin = M_ZERO 128 end if 129 130 ! load wavefunctions 131 call states_elec_load(restart, namespace, st, gr, ierr) 132 if(ierr /= 0) then 133 message(1) = "Unable to read wavefunctions." 134 call messages_fatal(1, namespace=namespace) 135 end if 136 137 POP_SUB(states_elec_look_and_load) 138 end subroutine states_elec_look_and_load 139 140 141 ! --------------------------------------------------------- 142 subroutine states_elec_dump(restart, st, gr, ierr, iter, lr, st_start_writing, verbose) 143 type(restart_t), intent(in) :: restart 144 type(states_elec_t), intent(in) :: st 145 type(grid_t), intent(in) :: gr 146 integer, intent(out) :: ierr 147 integer, optional, intent(in) :: iter 148 !> if this next argument is present, the lr wfs are stored instead of the gs wfs 149 type(lr_t), optional, intent(in) :: lr 150 integer, optional, intent(in) :: st_start_writing 151 logical, optional, intent(in) :: verbose 152 153 integer :: iunit_wfns, iunit_occs, iunit_states 154 integer :: err, err2(2), ik, idir, ist, idim, itot 155 integer :: root(1:P_STRATEGY_MAX) 156 character(len=MAX_PATH_LEN) :: filename 157 character(len=300) :: lines(3) 158 logical :: lr_wfns_are_associated, should_write, verbose_ 159 FLOAT :: kpoint(1:MAX_DIM) 160 FLOAT, allocatable :: dpsi(:), rff_global(:) 161 CMPLX, allocatable :: zpsi(:), zff_global(:) 162 163 PUSH_SUB(states_elec_dump) 164 165 verbose_ = optional_default(verbose, .true.) 166 167 ierr = 0 168 169 if(restart_skip(restart)) then 170 POP_SUB(states_elec_dump) 171 return 172 end if 173 174 if(verbose_) then 175 message(1) = "Info: Writing states." 176 call print_date(trim(message(1))//' ') 177 end if 178 179 call profiling_in(prof_write, "RESTART_WRITE") 180 181 if (present(lr)) then 182 lr_wfns_are_associated = (associated(lr%ddl_psi) .and. states_are_real(st)) .or. & 183 (associated(lr%zdl_psi) .and. states_are_complex(st)) 184 ASSERT(lr_wfns_are_associated) 185 end if 186 187 call restart_block_signals() 188 189 190 iunit_states = restart_open(restart, 'states') 191 write(lines(1), '(a20,1i10)') 'nst= ', st%nst 192 write(lines(2), '(a20,1i10)') 'dim= ', st%d%dim 193 write(lines(3), '(a20,1i10)') 'nik= ', st%d%nik 194 call restart_write(restart, iunit_states, lines, 3, err) 195 if (err /= 0) ierr = ierr + 1 196 call restart_close(restart, iunit_states) 197 198 199 iunit_wfns = restart_open(restart, 'wfns') 200 lines(1) = '# #k-point #st #dim filename' 201 if (states_are_real(st)) then 202 lines(2) = '%Real_Wavefunctions' 203 else 204 lines(2) = '%Complex_Wavefunctions' 205 end if 206 call restart_write(restart, iunit_wfns, lines, 2, err) 207 if (err /= 0) ierr = ierr + 2 208 209 210 iunit_occs = restart_open(restart, 'occs') 211 lines(1) = '# occupations | eigenvalue[a.u.] | Im(eigenvalue) [a.u.] | k-points | k-weights | filename | ik | ist | idim' 212 lines(2) = '%Occupations_Eigenvalues_K-Points' 213 call restart_write(restart, iunit_occs, lines, 2, err) 214 if (err /= 0) ierr = ierr + 4 215 216 217 if(states_are_real(st)) then 218 SAFE_ALLOCATE(dpsi(1:gr%mesh%np)) 219 SAFE_ALLOCATE(rff_global(1:gr%mesh%np_global)) 220 else 221 SAFE_ALLOCATE(zpsi(1:gr%mesh%np)) 222 SAFE_ALLOCATE(zff_global(1:gr%mesh%np_global)) 223 end if 224 225 itot = 1 226 root = -1 227 err2 = 0 228 do ik = 1, st%d%nik 229 kpoint = M_ZERO 230 kpoint(1:gr%sb%dim) = & 231 kpoints_get_point(gr%sb%kpoints, states_elec_dim_get_kpoint_index(st%d, ik), absolute_coordinates = .true.) 232 233 do ist = 1, st%nst 234 do idim = 1, st%d%dim 235 root(P_STRATEGY_DOMAINS) = mod(itot - 1, gr%mesh%mpi_grp%size) 236 write(filename,'(i10.10)') itot 237 238 write(lines(1), '(i8,a,i8,a,i8,3a)') ik, ' | ', ist, ' | ', idim, ' | "', trim(filename), '"' 239 call restart_write(restart, iunit_wfns, lines, 1, err) 240 if (err /= 0) err2(1) = err2(1) + 1 241 242 write(lines(1), '(e21.14,a,e21.14)') st%occ(ist,ik), ' | ', st%eigenval(ist, ik) 243 write(lines(1), '(a,a,e21.14)') trim(lines(1)), ' | ', CNST(0.0) 244 do idir = 1, gr%sb%dim 245 write(lines(1), '(a,a,e21.14)') trim(lines(1)), ' | ', kpoint(idir) 246 end do 247 write(lines(1), '(a,a,e21.14,a,i10.10,3(a,i8))') trim(lines(1)), & 248 ' | ', st%d%kweights(ik), ' | ', itot, ' | ', ik, ' | ', ist, ' | ', idim 249 call restart_write(restart, iunit_occs, lines, 1, err) 250 if (err /= 0) err2(1) = err2(1) + 1 251 252 should_write = st%st_start <= ist .and. ist <= st%st_end 253 if (should_write .and. present(st_start_writing)) then 254 if (ist < st_start_writing) should_write = .false. 255 end if 256 257 if (should_write) then 258 if (.not. present(lr)) then 259 if(st%d%kpt%start <= ik .and. ik <= st%d%kpt%end) then 260 if (states_are_real(st)) then 261 call states_elec_get_state(st, gr%mesh, idim, ist, ik, dpsi) 262 call drestart_write_mesh_function(restart, filename, gr%mesh, dpsi, err, root = root) 263 else 264 call states_elec_get_state(st, gr%mesh, idim, ist, ik, zpsi) 265 call zrestart_write_mesh_function(restart, filename, gr%mesh, zpsi, err, root = root) 266 end if 267 else 268 err = 0 269 end if 270 else 271 if(st%d%kpt%start <= ik .and. ik <= st%d%kpt%end) then 272 if (states_are_real(st)) then 273 call drestart_write_mesh_function(restart, filename, gr%mesh, & 274 lr%ddl_psi(:, idim, ist, ik), err, root = root) 275 else 276 call zrestart_write_mesh_function(restart, filename, gr%mesh, & 277 lr%zdl_psi(:, idim, ist, ik), err, root = root) 278 end if 279 else 280 err = 0 281 end if 282 end if 283 if (err /= 0) err2(2) = err2(2) + 1 284 end if 285 286 itot = itot + 1 287 end do ! st%d%dim 288 end do ! st%nst 289 end do ! st%d%nik 290 if (err2(1) /= 0) ierr = ierr + 8 291 if (err2(2) /= 0) ierr = ierr + 16 292 293 SAFE_DEALLOCATE_A(dpsi) 294 SAFE_DEALLOCATE_A(zpsi) 295 SAFE_DEALLOCATE_A(rff_global) 296 SAFE_DEALLOCATE_A(zff_global) 297 298 lines(1) = '%' 299 call restart_write(restart, iunit_occs, lines, 1, err) 300 if (err /= 0) ierr = ierr + 32 301 call restart_write(restart, iunit_wfns, lines, 1, err) 302 if (err /= 0) ierr = ierr + 64 303 if (present(iter)) then 304 write(lines(1),'(a,i7)') 'Iter = ', iter 305 call restart_write(restart, iunit_wfns, lines, 1, err) 306 if (err /= 0) ierr = ierr + 128 307 end if 308 309 call restart_close(restart, iunit_wfns) 310 call restart_close(restart, iunit_occs) 311 312 if(verbose_) then 313 message(1) = "Info: Finished writing states." 314 call print_date(trim(message(1))//' ') 315 end if 316 317 call restart_unblock_signals() 318 319 call profiling_out(prof_write) 320 POP_SUB(states_elec_dump) 321 return 322 end subroutine states_elec_dump 323 324 325 ! --------------------------------------------------------- 326 !> returns in ierr: 327 !! <0 => Fatal error, or nothing read 328 !! =0 => read all wavefunctions 329 !! >0 => could only read ierr wavefunctions 330 subroutine states_elec_load(restart, namespace, st, gr, ierr, iter, lr, lowest_missing, label, verbose, skip) 331 type(restart_t), intent(in) :: restart 332 type(namespace_t), intent(in) :: namespace 333 type(states_elec_t), intent(inout) :: st 334 type(grid_t), intent(in) :: gr 335 integer, intent(out) :: ierr 336 integer, optional, intent(out) :: iter 337 type(lr_t), optional, intent(inout) :: lr !< if present, the lr wfs are read instead of the gs wfs 338 integer, optional, intent(out) :: lowest_missing(:, :) !< all states below this one were read successfully 339 character(len=*), optional, intent(in) :: label 340 logical, optional, intent(in) :: verbose 341 logical, optional, intent(in) :: skip(:) !< A mask for reading or skipping the states. For expert use only. 342 343 integer :: states_elec_file, wfns_file, occ_file, err, ik, ist, idir, idim 344 integer :: idone, iread, ntodo 345 character(len=12) :: filename 346 character(len=1) :: char 347 logical, allocatable :: filled(:, :, :) 348 character(len=256) :: lines(3), label_ 349 character(len=50) :: str 350 351 FLOAT :: my_occ, imev, my_kweight 352 logical :: read_occ, lr_allocated, verbose_ 353 logical :: integral_occs 354 FLOAT, allocatable :: dpsi(:) 355 CMPLX, allocatable :: zpsi(:), zpsiL(:) 356 character(len=256), allocatable :: restart_file(:, :, :) 357 logical, allocatable :: restart_file_present(:, :, :) 358 FLOAT :: kpoint(MAX_DIM), read_kpoint(MAX_DIM) 359 360#if defined(HAVE_MPI) 361 integer :: iread_tmp 362 integer, allocatable :: lowest_missing_tmp(:, :) 363#endif 364 365 PUSH_SUB(states_elec_load) 366 367 ierr = 0 368 369 ! make sure these intent(out)`s are initialized no matter what 370 if (present(lowest_missing)) lowest_missing = 1 371 if (present(iter)) iter = 0 372 373 if (present(skip)) then 374 ASSERT(ubound(skip, dim = 1) == st%nst) 375 end if 376 377 if (restart_skip(restart)) then 378 ierr = -1 379 POP_SUB(states_elec_load) 380 return 381 end if 382 383 call profiling_in(prof_read, "RESTART_READ") 384 385 verbose_ = optional_default(verbose, .true.) 386 387 if (present(label)) then 388 label_ = trim(label) 389 else 390 if (present(lr)) then 391 label_ = " for linear response" 392 else 393 label_ = "" 394 end if 395 end if 396 397 message(1) = 'Info: Reading states' 398 if (len(trim(label_)) > 0) then 399 message(1) = trim(message(1)) // trim(label_) 400 end if 401 message(1) = trim(message(1)) // "." 402 if(verbose_) call print_date(trim(message(1))//' ') 403 404 if(.not. present(lr)) then 405 st%fromScratch = .false. ! obviously, we are using restart info 406 end if 407 408 ! If one restarts a GS calculation changing the %Occupations block, one 409 ! cannot read the occupations, otherwise these overwrite the ones from 410 ! the input file. restart_fixed_occ makes that we do use the ones in the file. 411 integral_occs = .true. ! only used if restart_fixed_occ 412 if (st%restart_fixed_occ) then 413 read_occ = .true. 414 st%fixed_occ = .true. 415 else 416 read_occ = .not. st%fixed_occ 417 end if 418 419 if(.not. present(lr)) then 420 st%eigenval(:, :) = M_ZERO 421 ! to be filled in from reading afterward 422 end if 423 424 if (.not. present(lr) .and. read_occ) then 425 st%occ(:, :) = M_ZERO 426 ! to be filled in from reading afterward 427 end if 428 429 ! sanity check 430 if (present(lr)) then 431 lr_allocated = (associated(lr%ddl_psi) .and. states_are_real(st)) .or. & 432 (associated(lr%zdl_psi) .and. states_are_complex(st)) 433 ASSERT(lr_allocated) 434 end if 435 436 states_elec_file = restart_open(restart, 'states') 437 ! sanity check on spin/k-points. Example file 'states': 438 ! nst= 2 439 ! dim= 1 440 ! nik= 2 441 call restart_read(restart, states_elec_file, lines, 3, err) 442 if (err /= 0) then 443 ierr = ierr - 2 444 else 445 read(lines(2), *) str, idim 446 read(lines(3), *) str, ik 447 if(idim == 2 .and. st%d%dim == 1) then 448 write(message(1),'(a)') 'Incompatible restart information: saved calculation is spinors, this one is not.' 449 call messages_warning(1, namespace=namespace) 450 ierr = ierr - 2**2 451 end if 452 if(idim == 1 .and. st%d%dim == 2) then 453 write(message(1),'(a)') 'Incompatible restart information: this calculation is spinors, saved one is not.' 454 call messages_warning(1, namespace=namespace) 455 ierr = ierr - 2**3 456 end if 457 if(ik < st%d%nik) then 458 write(message(1),'(a)') 'Incompatible restart information: not enough k-points.' 459 write(message(2),'(2(a,i6))') 'Expected ', st%d%nik, ' > Read ', ik 460 call messages_warning(2, namespace=namespace) 461 end if 462 ! We will check that they are the right k-points later, so we do not need to put a specific error here. 463 end if 464 call restart_close(restart, states_elec_file) 465 466 467 ! open files to read 468 wfns_file = restart_open(restart, 'wfns') 469 occ_file = restart_open(restart, 'occs') 470 call restart_read(restart, wfns_file, lines, 2, err) 471 if (err /= 0) then 472 ierr = ierr - 2**5 473 else if (states_are_real(st)) then 474 read(lines(2), '(a)') str 475 if (str(2:8) == 'Complex') then 476 message(1) = "Cannot read real states from complex wavefunctions." 477 call messages_warning(1, namespace=namespace) 478 ierr = ierr - 2**6 479 else if (str(2:5) /= 'Real') then 480 message(1) = "Restart file 'wfns' does not specify real/complex; cannot check compatibility." 481 call messages_warning(1, namespace=namespace) 482 end if 483 end if 484 ! complex can be restarted from real, so there is no problem. 485 486 ! Skip two lines. 487 call restart_read(restart, occ_file, lines, 2, err) 488 if (err /= 0) ierr = ierr - 2**7 489 490 ! If any error occured up to this point then it is not worth continuing, 491 ! as there something fundamentally wrong with the restart files 492 if (ierr /= 0) then 493 call restart_close(restart, wfns_file) 494 call restart_close(restart, occ_file) 495 call profiling_out(prof_read) 496 POP_SUB(states_elec_load) 497 return 498 end if 499 500 if (states_are_real(st)) then 501 SAFE_ALLOCATE(dpsi(1:gr%mesh%np)) 502 else 503 SAFE_ALLOCATE(zpsi(1:gr%mesh%np)) 504 end if 505 506 SAFE_ALLOCATE(restart_file(1:st%d%dim, st%st_start:st%st_end, 1:st%d%nik)) 507 SAFE_ALLOCATE(restart_file_present(1:st%d%dim, st%st_start:st%st_end, 1:st%d%nik)) 508 restart_file_present = .false. 509 510 ! Next we read the list of states from the files. 511 ! Errors in reading the information of a specific state from the files are ignored 512 ! at this point, because later we will skip reading the wavefunction of that state. 513 do 514 call restart_read(restart, wfns_file, lines, 1, err) 515 if (err == 0) then 516 read(lines(1), '(a)') char 517 if (char == '%') then 518 !We reached the end of the file 519 exit 520 else 521 read(lines(1), *) ik, char, ist, char, idim, char, filename 522 end if 523 end if 524 525 if (err /= 0 .or. index_is_wrong()) then 526 call restart_read(restart, occ_file, lines, 1, err) 527 cycle 528 end if 529 530 if (ist >= st%st_start .and. ist <= st%st_end .and. & 531 st%d%kpt%start <= ik .and. st%d%kpt%end >= ik) then 532 533 restart_file(idim, ist, ik) = trim(filename) 534 restart_file_present(idim, ist, ik) = .true. 535 end if 536 537 call restart_read(restart, occ_file, lines, 1, err) 538 if (.not. present(lr)) then ! do not read eigenvalues or occupations when reading linear response 539 ! # occupations | eigenvalue[a.u.] | Im(eigenvalue) [a.u.] | k-points | k-weights | filename | ik | ist | idim 540 541 if (err == 0) then 542 read(lines(1), *) my_occ, char, st%eigenval(ist, ik), char, imev, char, & 543 (read_kpoint(idir), char, idir = 1, gr%sb%dim), my_kweight 544 ! we do not want to read the k-weights, we have already set them appropriately 545 else 546 ! There is a problem with this states information, so we skip it. 547 restart_file_present(idim, ist, ik) = .false. 548 cycle 549 end if 550 551 kpoint(1:gr%sb%dim) = & 552 kpoints_get_point(gr%sb%kpoints, states_elec_dim_get_kpoint_index(st%d, ik), absolute_coordinates = .true.) 553 ! FIXME: maybe should ignore ik and just try to match actual vector k-points? 554 if (any(abs(kpoint(1:gr%sb%dim) - read_kpoint(1:gr%sb%dim)) > CNST(1e-12))) then 555 ! write only once for each k-point so as not to be too verbose 556 if (ist == 1) then 557 write(message(1),'(a,i6)') 'Incompatible restart information: k-point mismatch for ik ', ik 558 write(message(2),'(a,99f18.12)') ' Expected : ', kpoint(1:gr%sb%dim) 559 write(message(3),'(a,99f18.12)') ' Read : ', read_kpoint(1:gr%sb%dim) 560 call messages_warning(3, namespace=namespace) 561 end if 562 restart_file_present(idim, ist, ik) = .false. 563 end if 564 565 if (read_occ) then 566 st%occ(ist, ik) = my_occ 567 integral_occs = integral_occs .and. & 568 abs((st%occ(ist, ik) - st%smear%el_per_state) * st%occ(ist, ik)) <= M_EPSILON 569 end if 570 end if 571 end do 572 573 if (present(iter)) then 574 call restart_read(restart, wfns_file, lines, 1, err) 575 if (err /= 0) then 576 ierr = ierr - 2**8 577 else 578 read(lines(1), *) filename, filename, iter 579 end if 580 end if 581 582 call restart_close(restart, wfns_file) 583 call restart_close(restart, occ_file) 584 585 if (st%restart_fixed_occ) then 586 ! reset to overwrite whatever smearing may have been set earlier 587 call smear_init(st%smear, namespace, st%d%ispin, fixed_occ = .true., integral_occs = integral_occs, kpoints = gr%sb%kpoints) 588 end if 589 590 591 ! Now we read the wavefunctions. At this point we need to have all the information from the 592 ! states, occs, and wfns files in order to avoid serialisation of reads, as restart_read 593 ! works as a barrier. 594 595 SAFE_ALLOCATE(filled(1:st%d%dim, st%st_start:st%st_end, st%d%kpt%start:st%d%kpt%end)) 596 filled = .false. 597 598 if (present(lowest_missing)) lowest_missing = st%nst + 1 599 600 iread = 0 601 if (mpi_grp_is_root(mpi_world) .and. verbose_) then 602 idone = 0 603 ntodo = st%lnst*st%d%kpt%nlocal*st%d%dim 604 call loct_progress_bar(-1, ntodo) 605 end if 606 607 do ik = st%d%kpt%start, st%d%kpt%end 608 do ist = st%st_start, st%st_end 609 if(present(skip)) then 610 if(skip(ist)) cycle 611 end if 612 613 do idim = 1, st%d%dim 614 615 if (.not. restart_file_present(idim, ist, ik)) then 616 if (present(lowest_missing)) & 617 lowest_missing(idim, ik) = min(lowest_missing(idim, ik), ist) 618 cycle 619 end if 620 621 if (states_are_real(st)) then 622 call drestart_read_mesh_function(restart, restart_file(idim, ist, ik), gr%mesh, dpsi, err) 623 else 624 call zrestart_read_mesh_function(restart, restart_file(idim, ist, ik), gr%mesh, zpsi, err) 625 end if 626 627 if(states_are_real(st)) then 628 if(.not. present(lr)) then 629 call states_elec_set_state(st, gr%mesh, idim, ist, ik, dpsi) 630 else 631 call lalg_copy(gr%mesh%np, dpsi, lr%ddl_psi(:, idim, ist, ik)) 632 end if 633 else 634 if(.not. present(lr)) then 635 call states_elec_set_state(st, gr%mesh, idim, ist, ik, zpsi) 636 else 637 call lalg_copy(gr%mesh%np, zpsi, lr%zdl_psi(:, idim, ist, ik)) 638 end if 639 end if 640 641 642 if (err == 0) then 643 filled(idim, ist, ik) = .true. 644 iread = iread + 1 645 else if (present(lowest_missing)) then 646 lowest_missing(idim, ik) = min(lowest_missing(idim, ik), ist) 647 end if 648 649 if (mpi_grp_is_root(mpi_world) .and. verbose_) then 650 idone = idone + 1 651 call loct_progress_bar(idone, ntodo) 652 end if 653 654 end do 655 end do 656 end do 657 658 SAFE_DEALLOCATE_A(dpsi) 659 SAFE_DEALLOCATE_A(zpsi) 660 SAFE_DEALLOCATE_A(zpsiL) 661 SAFE_DEALLOCATE_A(restart_file) 662 SAFE_DEALLOCATE_A(restart_file_present) 663 664 if(mpi_grp_is_root(mpi_world) .and. verbose_) then 665 call messages_new_line() 666 end if 667 668#if defined(HAVE_MPI) 669 if(st%parallel_in_states .or. st%d%kpt%parallel) then 670 iread_tmp = iread 671 call MPI_Allreduce(iread_tmp, iread, 1, MPI_INTEGER, MPI_SUM, st%st_kpt_mpi_grp%comm, mpi_err) 672 end if 673 674 if(st%d%kpt%parallel) then 675 ! make sure all tasks know lowest_missing over all k-points 676 if(present(lowest_missing)) then 677 SAFE_ALLOCATE(lowest_missing_tmp(1:st%d%dim, 1:st%d%nik)) 678 lowest_missing_tmp = lowest_missing 679 call MPI_Allreduce(lowest_missing_tmp, lowest_missing, st%d%dim*st%d%nik, & 680 MPI_INTEGER, MPI_MIN, st%d%kpt%mpi_grp%comm, mpi_err) 681 SAFE_DEALLOCATE_A(lowest_missing_tmp) 682 end if 683 end if 684#endif 685 686 if (.not. present(lr) .and. .not. present(skip)) call fill_random() 687 ! it is better to initialize lr wfns to zero 688 689 SAFE_DEALLOCATE_A(filled) 690 691 if (ierr == 0 .and. iread /= st%nst * st%d%nik * st%d%dim) then 692 if(iread > 0) then 693 ierr = iread 694 else 695 ierr = -1 696 end if 697 ! otherwise ierr = 0 would mean either all was read correctly, or nothing at all was read! 698 699 if(.not. present(lr)) then 700 write(str, '(a,i5)') 'Reading states.' 701 else 702 write(str, '(a,i5)') 'Reading states information for linear response.' 703 end if 704 call messages_print_stress(stdout, trim(str), namespace=namespace) 705 if(.not. present(skip)) then 706 write(message(1),'(a,i6,a,i6,a)') 'Only ', iread,' files out of ', & 707 st%nst * st%d%nik * st%d%dim, ' could be read.' 708 else 709 write(message(1),'(a,i6,a,i6,a)') 'Only ', iread,' files out of ', & 710 st%nst * st%d%nik * st%d%dim, ' were loaded.' 711 ierr = 0 712 end if 713 call messages_info(1) 714 call messages_print_stress(stdout, namespace=namespace) 715 end if 716 717 message(1) = 'Info: States reading done.' 718 if(verbose_) call print_date(trim(message(1))//' ') 719 720 call profiling_out(prof_read) 721 POP_SUB(states_elec_load) 722 723 contains 724 725 ! --------------------------------------------------------- 726 subroutine fill_random() !< Put random function in orbitals that could not be read. 727 PUSH_SUB(states_elec_load.fill_random) 728 729 do ik = st%d%kpt%start, st%d%kpt%end 730 731 do ist = st%st_start, st%st_end 732 do idim = 1, st%d%dim 733 if(filled(idim, ist, ik)) cycle 734 735 call states_elec_generate_random(st, gr%mesh, gr%sb, ist, ist, ik, ik) 736 end do 737 end do 738 end do 739 740 POP_SUB(states_elec_load.fill_random) 741 end subroutine fill_random 742 743 ! --------------------------------------------------------- 744 745 logical function index_is_wrong() !< .true. if the index (idim, ist, ik) is not present in st structure... 746 PUSH_SUB(states_elec_load.index_is_wrong) 747 748 if(idim > st%d%dim .or. idim < 1 .or. & 749 ist > st%nst .or. ist < 1 .or. & 750 ik > st%d%nik .or. ik < 1) then 751 index_is_wrong = .true. 752 else 753 index_is_wrong = .false. 754 end if 755 756 POP_SUB(states_elec_load.index_is_wrong) 757 end function index_is_wrong 758 759 end subroutine states_elec_load 760 761 762 subroutine states_elec_dump_rho(restart, st, gr, ierr, iter) 763 type(restart_t), intent(in) :: restart 764 type(states_elec_t), intent(in) :: st 765 type(grid_t), intent(in) :: gr 766 integer, intent(out) :: ierr 767 integer, optional, intent(in) :: iter 768 769 integer :: iunit, isp, err, err2(2) 770 character(len=80) :: filename 771 character(len=300) :: lines(2) 772 FLOAT, pointer :: rho(:), rho_fine(:) 773 774 PUSH_SUB(states_elec_dump_rho) 775 776 ierr = 0 777 778 if (restart_skip(restart)) then 779 POP_SUB(states_elec_dump_rho) 780 return 781 end if 782 783 if (debug%info) then 784 message(1) = "Debug: Writing density restart." 785 call messages_info(1) 786 end if 787 788 call restart_block_signals() 789 790 !write the densities 791 iunit = restart_open(restart, 'density') 792 lines(1) = '# #spin #nspin filename' 793 lines(2) = '%densities' 794 call restart_write(restart, iunit, lines, 2, err) 795 if (err /= 0) ierr = ierr + 1 796 797 if(gr%have_fine_mesh) then 798 SAFE_ALLOCATE(rho(1:gr%mesh%np)) 799 SAFE_ALLOCATE(rho_fine(1:gr%fine%mesh%np)) 800 end if 801 802 err2 = 0 803 do isp = 1, st%d%nspin 804 if(st%d%nspin==1) then 805 write(filename, fmt='(a)') 'density' 806 else 807 write(filename, fmt='(a,i1)') 'density-sp', isp 808 end if 809 write(lines(1), '(i8,a,i8,a)') isp, ' | ', st%d%nspin, ' | "'//trim(adjustl(filename))//'"' 810 call restart_write(restart, iunit, lines, 1, err) 811 if (err /= 0) err2(1) = err2(1) + 1 812 813 if(gr%have_fine_mesh)then 814 rho_fine(1:gr%fine%mesh%np) = st%rho(1:gr%fine%mesh%np,isp) 815 call dmultigrid_fine2coarse(gr%fine%tt, gr%fine%der, gr%mesh, rho_fine, rho, INJECTION) 816 call drestart_write_mesh_function(restart, filename, gr%mesh, rho, err) 817 else 818 call drestart_write_mesh_function(restart, filename, gr%mesh, st%rho(:,isp), err) 819 end if 820 if (err /= 0) err2(2) = err2(2) + 1 821 822 end do 823 if (err2(1) /= 0) ierr = ierr + 2 824 if (err2(2) /= 0) ierr = ierr + 4 825 826 if(gr%have_fine_mesh)then 827 SAFE_DEALLOCATE_P(rho) 828 SAFE_DEALLOCATE_P(rho_fine) 829 end if 830 831 lines(1) = '%' 832 call restart_write(restart, iunit, lines, 1, err) 833 if (err /= 0) ierr = ierr + 8 834 if (present(iter)) then 835 write(lines(1),'(a,i7)') 'Iter = ', iter 836 call restart_write(restart, iunit, lines, 1, err) 837 if (err /= 0) ierr = ierr + 16 838 end if 839 call restart_close(restart, iunit) 840 841 call restart_unblock_signals() 842 843 if (debug%info) then 844 message(1) = "Debug: Writing density restart done." 845 call messages_info(1) 846 end if 847 848 POP_SUB(states_elec_dump_rho) 849 end subroutine states_elec_dump_rho 850 851 852 ! --------------------------------------------------------- 853 subroutine states_elec_load_rho(restart, st, gr, ierr) 854 type(restart_t), intent(in) :: restart 855 type(states_elec_t), intent(inout) :: st 856 type(grid_t), intent(in) :: gr 857 integer, intent(out) :: ierr 858 859 integer :: err, err2, isp 860 character(len=12) :: filename 861 FLOAT, allocatable :: rho_coarse(:) 862 863 PUSH_SUB(states_elec_load_rho) 864 865 ierr = 0 866 867 if (restart_skip(restart)) then 868 ierr = -1 869 POP_SUB(states_elec_load_rho) 870 return 871 end if 872 873 if (debug%info) then 874 message(1) = "Debug: Reading density restart." 875 call messages_info(1) 876 end if 877 878 ! skip for now, since we know what the files are going to be called 879 !read the densities 880! iunit_rho = io_open(trim(dir)//'/density', restart%namespace, action='write') 881! call iopar_read(st%dom_st_kpt_mpi_grp, iunit_rho, line, err) 882! call iopar_read(st%dom_st_kpt_mpi_grp, iunit_rho, line, err) 883! we could read the iteration 'iter' too, not sure if that is useful. 884 885 if(gr%have_fine_mesh) then 886 SAFE_ALLOCATE(rho_coarse(1:gr%mesh%np_part)) 887 end if 888 889 err2 = 0 890 do isp = 1, st%d%nspin 891 if(st%d%nspin==1) then 892 write(filename, fmt='(a)') 'density' 893 else 894 write(filename, fmt='(a,i1)') 'density-sp', isp 895 end if 896! if(mpi_grp_is_root(st%dom_st_kpt_mpi_grp)) then 897! read(iunit_rho, '(i8,a,i8,a)') isp, ' | ', st%d%nspin, ' | "'//trim(adjustl(filename))//'"' 898! end if 899 if(gr%have_fine_mesh)then 900 call drestart_read_mesh_function(restart, filename, gr%mesh, rho_coarse, err) 901 call dmultigrid_coarse2fine(gr%fine%tt, gr%der, gr%fine%mesh, rho_coarse, st%rho(:,isp), order = 2) 902 else 903 call drestart_read_mesh_function(restart, filename, gr%mesh, st%rho(:,isp), err) 904 end if 905 if (err /= 0) err2 = err2 + 1 906 907 end do 908 if (err2 /= 0) ierr = ierr + 1 909 910 if(gr%have_fine_mesh)then 911 SAFE_DEALLOCATE_A(rho_coarse) 912 end if 913 914 if (debug%info) then 915 message(1) = "Debug: Reading density restart done." 916 call messages_info(1) 917 end if 918 919 POP_SUB(states_elec_load_rho) 920 end subroutine states_elec_load_rho 921 922 subroutine states_elec_dump_frozen(restart, st, gr, ierr) 923 type(restart_t), intent(in) :: restart 924 type(states_elec_t), intent(in) :: st 925 type(grid_t), intent(in) :: gr 926 integer, intent(out) :: ierr 927 928 integer :: isp, err, err2(2), idir 929 character(len=80) :: filename 930 931 PUSH_SUB(states_elec_dump_frozen) 932 933 ierr = 0 934 935 ASSERT(associated(st%frozen_rho)) 936 937 if (restart_skip(restart)) then 938 POP_SUB(states_elec_dump_frozen) 939 return 940 end if 941 942 if (debug%info) then 943 message(1) = "Debug: Writing frozen densities restart." 944 call messages_info(1) 945 end if 946 947 call restart_block_signals() 948 949 err2 = 0 950 do isp = 1, st%d%nspin 951 if(st%d%nspin==1) then 952 write(filename, fmt='(a)') 'frozen_rho' 953 else 954 write(filename, fmt='(a,i1)') 'frozen_rho-sp', isp 955 end if 956 957 call drestart_write_mesh_function(restart, filename, gr%mesh, st%frozen_rho(:,isp), err) 958 if (err /= 0) err2(2) = err2(2) + 1 959 960 if(associated(st%frozen_tau)) then 961 if(st%d%nspin==1) then 962 write(filename, fmt='(a)') 'frozen_tau' 963 else 964 write(filename, fmt='(a,i1)') 'frozen_tau-sp', isp 965 end if 966 call drestart_write_mesh_function(restart, filename, gr%mesh, st%frozen_tau(:,isp), err) 967 if (err /= 0) err2 = err2 + 1 968 end if 969 970 if(associated(st%frozen_gdens)) then 971 do idir = 1, gr%sb%dim 972 if(st%d%nspin==1) then 973 write(filename, fmt='(a,i1)') 'frozen_gdens-dir', idir 974 else 975 write(filename, fmt='(a,i1,a,i1)') 'frozen_tau-dir', idir, '-', isp 976 end if 977 call drestart_write_mesh_function(restart, filename, gr%mesh, st%frozen_gdens(:,idir,isp), err) 978 if (err /= 0) err2 = err2 + 1 979 end do 980 end if 981 982 if(associated(st%frozen_ldens)) then 983 if(st%d%nspin==1) then 984 write(filename, fmt='(a)') 'frozen_ldens' 985 else 986 write(filename, fmt='(a,i1)') 'frozen_ldens-sp', isp 987 end if 988 call drestart_write_mesh_function(restart, filename, gr%mesh, st%frozen_ldens(:,isp), err) 989 if (err /= 0) err2 = err2 + 1 990 end if 991 992 end do 993 if (err2(1) /= 0) ierr = ierr + 2 994 if (err2(2) /= 0) ierr = ierr + 4 995 996 call restart_unblock_signals() 997 998 if (debug%info) then 999 message(1) = "Debug: Writing frozen densities restart done." 1000 call messages_info(1) 1001 end if 1002 1003 POP_SUB(states_elec_dump_frozen) 1004 end subroutine states_elec_dump_frozen 1005 1006 1007 ! --------------------------------------------------------- 1008 subroutine states_elec_load_frozen(restart, st, gr, ierr) 1009 type(restart_t), intent(in) :: restart 1010 type(states_elec_t), intent(inout) :: st 1011 type(grid_t), intent(in) :: gr 1012 integer, intent(out) :: ierr 1013 1014 integer :: err, err2, isp, idir 1015 character(len=80) :: filename 1016 1017 PUSH_SUB(states_elec_load_frozen) 1018 1019 ASSERT(associated(st%frozen_rho)) 1020 1021 ierr = 0 1022 1023 if (restart_skip(restart)) then 1024 ierr = -1 1025 POP_SUB(states_elec_load_frozen) 1026 return 1027 end if 1028 1029 if (debug%info) then 1030 message(1) = "Debug: Reading densities restart." 1031 call messages_info(1) 1032 end if 1033 1034 err2 = 0 1035 do isp = 1, st%d%nspin 1036 if(st%d%nspin==1) then 1037 write(filename, fmt='(a)') 'frozen_rho' 1038 else 1039 write(filename, fmt='(a,i1)') 'frozen_rho-sp', isp 1040 end if 1041 call drestart_read_mesh_function(restart, filename, gr%mesh, st%frozen_rho(:,isp), err) 1042 if (err /= 0) err2 = err2 + 1 1043 1044 if(associated(st%frozen_tau)) then 1045 if(st%d%nspin==1) then 1046 write(filename, fmt='(a)') 'frozen_tau' 1047 else 1048 write(filename, fmt='(a,i1)') 'frozen_tau-sp', isp 1049 end if 1050 call drestart_read_mesh_function(restart, filename, gr%mesh, st%frozen_tau(:,isp), err) 1051 if (err /= 0) err2 = err2 + 1 1052 end if 1053 1054 if(associated(st%frozen_gdens)) then 1055 do idir = 1, gr%sb%dim 1056 if(st%d%nspin==1) then 1057 write(filename, fmt='(a,i1)') 'frozen_gdens-dir', idir 1058 else 1059 write(filename, fmt='(a,i1,a,i1)') 'frozen_tau-dir', idir, '-', isp 1060 end if 1061 call drestart_read_mesh_function(restart, filename, gr%mesh, st%frozen_gdens(:,idir,isp), err) 1062 if (err /= 0) err2 = err2 + 1 1063 end do 1064 end if 1065 1066 if(associated(st%frozen_ldens)) then 1067 if(st%d%nspin==1) then 1068 write(filename, fmt='(a)') 'frozen_ldens' 1069 else 1070 write(filename, fmt='(a,i1)') 'frozen_ldens-sp', isp 1071 end if 1072 call drestart_read_mesh_function(restart, filename, gr%mesh, st%frozen_ldens(:,isp), err) 1073 if (err /= 0) err2 = err2 + 1 1074 end if 1075 1076 end do 1077 if (err2 /= 0) ierr = ierr + 1 1078 1079 if (debug%info) then 1080 message(1) = "Debug: Reading frozen densities restart done." 1081 call messages_info(1) 1082 end if 1083 1084 POP_SUB(states_elec_load_frozen) 1085 end subroutine states_elec_load_frozen 1086 1087 1088 ! --------------------------------------------------------- 1089 !> the routine reads formulas for user-defined wavefunctions 1090 !! from the input file and fills the respective orbitals 1091 subroutine states_elec_read_user_def_orbitals(mesh, namespace, st) 1092 type(mesh_t), intent(in) :: mesh 1093 type(namespace_t), intent(in) :: namespace 1094 type(states_elec_t), intent(inout) :: st 1095 1096 type(block_t) :: blk 1097 integer :: ip, id, is, ik, nstates, state_from, ierr, ncols 1098 integer :: ib, idim, inst, inik, normalize 1099 FLOAT :: xx(MAX_DIM), rr, psi_re, psi_im 1100 character(len=150) :: filename 1101 CMPLX, allocatable :: zpsi(:, :) 1102 1103 integer, parameter :: & 1104 STATE_FROM_FORMULA = 1, & 1105 STATE_FROM_FILE = -10010, & 1106 NORMALIZE_YES = 1, & 1107 NORMALIZE_NO = 0 1108 1109 PUSH_SUB(states_elec_read_user_def_orbitals) 1110 1111 !%Variable UserDefinedStates 1112 !%Type block 1113 !%Section States 1114 !%Description 1115 !% Instead of using the ground state as initial state for 1116 !% time-propagations it might be interesting in some cases 1117 !% to specify alternate states. Like with user-defined 1118 !% potentials, this block allows you to specify formulas for 1119 !% the orbitals at <i>t</i>=0. 1120 !% 1121 !% Example: 1122 !% 1123 !% <tt>%UserDefinedStates 1124 !% <br> 1 | 1 | 1 | formula | "exp(-r^2)*exp(-i*0.2*x)" | normalize_yes 1125 !% <br>%</tt> 1126 !% 1127 !% The first column specifies the component of the spinor, 1128 !% the second column the number of the state and the third 1129 !% contains <i>k</i>-point and spin quantum numbers. Column four 1130 !% indicates that column five should be interpreted as a formula 1131 !% for the corresponding orbital. 1132 !% 1133 !% Alternatively, if column four states <tt>file</tt> the state will 1134 !% be read from the file given in column five. 1135 !% 1136 !% <tt>%UserDefinedStates 1137 !% <br> 1 | 1 | 1 | file | "/path/to/file" | normalize_no 1138 !% <br>%</tt> 1139 !% 1140 !% Octopus reads first the ground-state orbitals from 1141 !% the <tt>restart/gs</tt> directory. Only the states that are 1142 !% specified in the above block will be overwritten with 1143 !% the given analytic expression for the orbital. 1144 !% 1145 !% The sixth (optional) column indicates whether <tt>Octopus</tt> should renormalize 1146 !% the orbital. The default (no sixth column given) is to renormalize. 1147 !% 1148 !%Option file -10010 1149 !% Read initial orbital from file. 1150 !% Accepted file formats, detected by extension: obf, ncdf and csv (real only). 1151 !%Option formula 1 1152 !% Calculate initial orbital by given analytic expression. 1153 !%Option normalize_yes 1 1154 !% Normalize orbitals (default). 1155 !%Option normalize_no 0 1156 !% Do not normalize orbitals. 1157 !%End 1158 if(parse_block(namespace, 'UserDefinedStates', blk) == 0) then 1159 1160 call messages_print_stress(stdout, trim('Substitution of orbitals'), namespace=namespace) 1161 1162 ! find out how many lines (i.e. states) the block has 1163 nstates = parse_block_n(blk) 1164 1165 SAFE_ALLOCATE(zpsi(1:mesh%np, 1:st%d%dim)) 1166 1167 ! read all lines 1168 do ib = 1, nstates 1169 ! Check that number of columns is five or six. 1170 ncols = parse_block_cols(blk, ib - 1) 1171 if(ncols < 5 .or. ncols > 6) then 1172 message(1) = 'Each line in the UserDefinedStates block must have' 1173 message(2) = 'five or six columns.' 1174 call messages_fatal(2, namespace=namespace) 1175 end if 1176 1177 call parse_block_integer(blk, ib - 1, 0, idim) 1178 call parse_block_integer(blk, ib - 1, 1, inst) 1179 call parse_block_integer(blk, ib - 1, 2, inik) 1180 1181 ! Calculate from expression or read from file? 1182 call parse_block_integer(blk, ib - 1, 3, state_from) 1183 1184 ! loop over all states 1185 do id = 1, st%d%dim 1186 do is = 1, st%nst 1187 do ik = 1, st%d%nik 1188 1189 ! does the block entry match and is this node responsible? 1190 if(.not.(id == idim .and. is == inst .and. ik == inik & 1191 .and. st%st_start <= is .and. st%st_end >= is & 1192 .and. st%d%kpt%start <= ik .and. st%d%kpt%end >= ik) ) cycle 1193 1194 select case(state_from) 1195 1196 case(STATE_FROM_FORMULA) 1197 ! parse formula string 1198 call parse_block_string( & 1199 blk, ib - 1, 4, st%user_def_states(id, is, ik)) 1200 1201 write(message(1), '(a,3i5)') 'Substituting state of orbital with k, ist, dim = ', ik, is, id 1202 write(message(2), '(2a)') ' with the expression:' 1203 write(message(3), '(2a)') ' ',trim(st%user_def_states(id, is, ik)) 1204 call messages_info(3) 1205 1206 ! convert to C string 1207 call conv_to_C_string(st%user_def_states(id, is, ik)) 1208 1209 ! fill states with user-defined formulas 1210 do ip = 1, mesh%np 1211 xx = mesh%x(ip, :) 1212 rr = sqrt(sum(xx(:)**2)) 1213 1214 ! parse user-defined expressions 1215 call parse_expression(psi_re, psi_im, mesh%sb%dim, xx, rr, M_ZERO, st%user_def_states(id, is, ik)) 1216 ! fill state 1217 zpsi(ip, 1) = psi_re + M_zI * psi_im 1218 end do 1219 1220 case(STATE_FROM_FILE) 1221 ! The input format can be coded in column four now. As it is 1222 ! not used now, we just say "file". 1223 ! Read the filename. 1224 call parse_block_string(blk, ib - 1, 4, filename) 1225 1226 write(message(1), '(a,3i5)') 'Substituting state of orbital with k, ist, dim = ', ik, is, id 1227 write(message(2), '(2a)') ' with data from file:' 1228 write(message(3), '(2a)') ' ',trim(filename) 1229 call messages_info(3) 1230 1231 ! finally read the state 1232 call zio_function_input(filename, namespace, mesh, zpsi(:, 1), ierr) 1233 if (ierr > 0) then 1234 message(1) = 'Could not read the file!' 1235 write(message(2),'(a,i1)') 'Error code: ', ierr 1236 call messages_fatal(2, namespace=namespace) 1237 end if 1238 1239 case default 1240 message(1) = 'Wrong entry in UserDefinedStates, column 4.' 1241 message(2) = 'You may state "formula" or "file" here.' 1242 call messages_fatal(2, namespace=namespace) 1243 end select 1244 1245 call states_elec_set_state(st, mesh, id, is, ik, zpsi(:, 1)) 1246 1247 ! normalize orbital 1248 if(parse_block_cols(blk, ib - 1) == 6) then 1249 call parse_block_integer(blk, ib - 1, 5, normalize) 1250 else 1251 normalize = 1 1252 end if 1253 select case(normalize) 1254 case(NORMALIZE_NO) 1255 case(NORMALIZE_YES) 1256 call states_elec_get_state(st, mesh, is, ik, zpsi) 1257 call zmf_normalize(mesh, st%d%dim, zpsi) 1258 call states_elec_set_state(st, mesh, is, ik, zpsi) 1259 case default 1260 message(1) = 'The sixth column in UserDefinedStates may either be' 1261 message(2) = '"normalize_yes" or "normalize_no"' 1262 call messages_fatal(2, namespace=namespace) 1263 end select 1264 1265 end do 1266 end do 1267 end do 1268 1269 end do 1270 1271 SAFE_DEALLOCATE_A(zpsi) 1272 1273 call parse_block_end(blk) 1274 call messages_print_stress(stdout, namespace=namespace) 1275 1276 else 1277 message(1) = "'UserDefinedStates' has to be specified as block." 1278 call messages_fatal(1, namespace=namespace) 1279 end if 1280 1281 POP_SUB(states_elec_read_user_def_orbitals) 1282 end subroutine states_elec_read_user_def_orbitals 1283 1284 ! --------------------------------------------------------- 1285 ! This is needed as for the generalized Bloch theorem we need to label 1286 ! the states from the expectation value of Sz computed from the GS. 1287 subroutine states_elec_dump_spin(restart, st, ierr) 1288 type(restart_t), intent(in) :: restart 1289 type(states_elec_t), intent(in) :: st 1290 integer, intent(out) :: ierr 1291 1292 integer :: iunit_spin 1293 integer :: err, err2(2), ik, ist 1294 character(len=300) :: lines(3) 1295 1296 PUSH_SUB(states_elec_dump_spin) 1297 1298 ierr = 0 1299 1300 if(restart_skip(restart)) then 1301 POP_SUB(states_dump_spin) 1302 return 1303 end if 1304 1305 call profiling_in(prof_write, "RESTART_WRITE_SPIN") 1306 1307 call restart_block_signals() 1308 1309 iunit_spin = restart_open(restart, 'spin') 1310 lines(1) = '# #k-point #st #spin(x) spin(y) spin(z)' 1311 call restart_write(restart, iunit_spin, lines, 1, err) 1312 if (err /= 0) ierr = ierr + 1 1313 1314 err2 = 0 1315 do ik = 1, st%d%nik 1316 do ist = 1, st%nst 1317 write(lines(1), '(i8,a,i8,3(a,f18.12))') ik, ' | ', ist, ' | ', & 1318 st%spin(1,ist,ik), ' | ', st%spin(2,ist,ik),' | ', st%spin(3,ist,ik) 1319 call restart_write(restart, iunit_spin, lines, 1, err) 1320 if (err /= 0) err2(1) = err2(1) + 1 1321 end do ! st%nst 1322 end do ! st%d%nik 1323 lines(1) = '%' 1324 call restart_write(restart, iunit_spin, lines, 1, err) 1325 if (err2(1) /= 0) ierr = ierr + 8 1326 if (err2(2) /= 0) ierr = ierr + 16 1327 1328 call restart_close(restart, iunit_spin) 1329 1330 call restart_unblock_signals() 1331 1332 call profiling_out(prof_write) 1333 POP_SUB(states_elec_dump_spin) 1334 end subroutine states_elec_dump_spin 1335 1336 1337 1338 1339 ! --------------------------------------------------------- 1340 !> returns in ierr: 1341 !! <0 => Fatal error, or nothing read 1342 !! =0 => read all wavefunctions 1343 !! >0 => could only read ierr wavefunctions 1344 subroutine states_elec_load_spin(restart, st, ierr) 1345 type(restart_t), intent(in) :: restart 1346 type(states_elec_t), intent(inout) :: st 1347 integer, intent(out) :: ierr 1348 1349 integer :: spin_file, err, ik, ist 1350 character(len=256) :: lines(3) 1351 FLOAT :: spin(3) 1352 character(len=1) :: char 1353 1354 1355 PUSH_SUB(states_elec_load_spin) 1356 1357 ierr = 0 1358 1359 if (restart_skip(restart)) then 1360 ierr = -1 1361 POP_SUB(states_load_spin) 1362 return 1363 end if 1364 1365 call profiling_in(prof_read, "RESTART_READ_SPIN") 1366 1367 ! open files to read 1368 spin_file = restart_open(restart, 'spin') 1369 ! Skip two lines. 1370 call restart_read(restart, spin_file, lines, 1, err) 1371 if (err /= 0) ierr = ierr - 2**7 1372 1373 ! If any error occured up to this point then it is not worth continuing, 1374 ! as there something fundamentally wrong with the restart files 1375 if (ierr /= 0) then 1376 call restart_close(restart, spin_file) 1377 call profiling_out(prof_read) 1378 POP_SUB(states_load_spin) 1379 return 1380 end if 1381 1382 ! Next we read the list of states from the files. 1383 ! Errors in reading the information of a specific state from the files are ignored 1384 ! at this point, because later we will skip reading the wavefunction of that state. 1385 do 1386 call restart_read(restart, spin_file, lines, 1, err) 1387 read(lines(1), '(a)') char 1388 if (char == '%') then 1389 !We reached the end of the file 1390 exit 1391 else 1392 read(lines(1), *) ik, char, ist, char, spin(1), char, spin(2), char, spin(3) 1393 end if 1394 1395 if (err /= 0) cycle 1396 1397 st%spin(1:3, ist, ik) = spin(1:3) 1398 end do 1399 1400 call restart_close(restart, spin_file) 1401 1402 call profiling_out(prof_read) 1403 POP_SUB(states_elec_load_spin) 1404 end subroutine states_elec_load_spin 1405 1406end module states_elec_restart_oct_m 1407 1408 1409!! Local Variables: 1410!! mode: f90 1411!! coding: utf-8 1412!! End: 1413