1!! Copyright (C) 2019 R. Jestaedt, F. Bonafe, H. Appel, A. Rubio 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_mxll_restart_oct_m 22 use global_oct_m 23 use grid_oct_m 24 use io_function_oct_m 25 use lalg_basic_oct_m 26 use loct_oct_m 27 use mesh_oct_m 28 use mesh_function_oct_m 29 use messages_oct_m 30 use mpi_oct_m 31 use multicomm_oct_m 32 use multigrid_oct_m 33 use namespace_oct_m 34 use parser_oct_m 35 use profiling_oct_m 36 use restart_oct_m 37 use simul_box_oct_m 38 use states_mxll_oct_m 39 use states_elec_dim_oct_m 40 use string_oct_m 41 use types_oct_m 42 use unit_system_oct_m 43 use unit_oct_m 44 use states_elec_restart_oct_m 45 use batch_oct_m 46 use batch_ops_oct_m 47 48 implicit none 49 50 type(profile_t), save :: prof_read, prof_write 51 52 private 53 public :: & 54 states_mxll_read_user_def, & 55 states_mxll_dump, & 56 states_mxll_load 57 58contains 59 60 !---------------------------------------------------------- 61 subroutine states_mxll_read_user_def(mesh, st, user_def_rs_state, namespace) 62 type(mesh_t), intent(inout) :: mesh 63 type(states_mxll_t), intent(inout) :: st 64 type(namespace_t), intent(in) :: namespace 65 CMPLX, intent(inout) :: user_def_rs_state(:,:) 66 67 type(block_t) :: blk 68 integer :: il, nlines, idim, ncols, ip, state_from, ierr, maxwell_field 69 FLOAT :: xx(MAX_DIM), rr, e_value, dummy, b_value 70 FLOAT, allocatable :: e_field(:), b_field(:) 71 CMPLX, allocatable :: rs_state(:,:), rs_state_add(:,:) 72 character(len=150), pointer :: filename_e_field, filename_b_field 73 character(1) :: cdim 74 75 integer, parameter :: & 76 STATE_FROM_FORMULA = 1, & 77 STATE_FROM_FILE = -10010 78 79 PUSH_SUB(states_mxll_read_user_def) 80 81 SAFE_ALLOCATE(rs_state(1:mesh%np_part,1:st%dim)) 82 SAFE_ALLOCATE(rs_state_add(1:mesh%np_part,1:st%dim)) 83 84 ! Set electromagnetic field equal to zero in the whole simulation box. 85 rs_state(:,:) = M_ZERO 86 rs_state_add(:,:) = M_ZERO 87 88 !%Variable UserDefinedInitialMaxwellStates 89 !%Type block 90 !%Section MaxwellStates 91 !%Description 92 !% The initial electromagnetic fields can be set by the user 93 !% with the <tt>UserDefinedMaxwellStates</tt> block variable. 94 !% The electromagnetic fields have to fulfill the 95 !% Maxwells equations in vacuum. 96 !% 97 !% Example: 98 !% 99 !% <tt>%UserDefinedMaxwellStates 100 !% <br> 2 | formula | "magnetic_field" | "-1/P_c * sin(x)" 101 !% <br> 3 | formula | "electric_field" | " sin(x) " 102 !% <br>%</tt> 103 !% 104 !% The first column specifies the component of the dimension of 105 !% the electric field and magnetic field. Column four 106 !% indicates that column five and six should be interpreted 107 !% as a formula for the corresponding orbital. P_c is the 108 !% speed of light constant. 109 !% 110 !% Alternatively, if column four states <tt>file</tt> the 111 !% electric field and magnetic field will be read from 112 !% the files given in column five and six. 113 !% 114 !% <tt>%MUserDefinedMaxwellStates 115 !% <br> 3 | file | electric_field | "/path/to/file_electric_field_of_dimension_3" 116 !% <br> 2 | file | magnetic_field | "/path/to/file_magnetic_field_of_dimension_2" 117 !% <br>%</tt> 118 !% 119 !%Option file -10010 120 !% Read initial orbital from file. 121 !% Accepted file formats: obf, ncdf and csv. 122 !%Option formula 1 123 !% Calculate initial orbital by given analytic expression. 124 !%Option electric_field 1 125 !% This row defines the electric field component of the corresponding dimension 126 !%Option magnetic_field 2 127 !% This row defines the magnetic field component of the corresponding dimension 128 !%End 129 130 if(parse_block(namespace, 'UserDefinedInitialMaxwellStates', blk) == 0) then 131 132 ! find out how many lines (i.e. states) the block has 133 nlines = parse_block_n(blk) 134 135 write(message(1), '(a,i5)') 'Maxwell electromagnetic fields are added.' 136 write(message(2), '(a,i5)') '' 137 call messages_info(2) 138 139 ! read all lines 140 do il = 1, nlines 141 ! Check that number of columns is five or six. 142 ncols = parse_block_cols(blk, il - 1) 143 if(ncols < 4 .or. ncols > 4) then 144 message(1) = 'Each line in the UserDefinedMaxwellStates block must have' 145 message(2) = 'four columns.' 146 call messages_fatal(2, namespace=namespace) 147 end if 148 149 call parse_block_integer(blk, il - 1, 0, idim) 150 write(cdim,'(I1)')idim 151 152 ! Calculate from expression or read from file? 153 call parse_block_integer(blk, il - 1, 1, state_from) 154 155 select case(state_from) 156 157 case(STATE_FROM_FORMULA) 158 ! parse formula string 159 call parse_block_integer(blk, il - 1, 2, maxwell_field ) 160 if (maxwell_field == OPTION__USERDEFINEDINITIALMAXWELLSTATES__ELECTRIC_FIELD) then 161 call parse_block_string( blk, il - 1, 3, st%user_def_e_field(idim)) 162 call messages_write(" E-field in dimension "//trim(cdim)//" : "//trim(st%user_def_e_field(idim)), fmt='(a,i1,2a)') 163 call conv_to_C_string(st%user_def_e_field(idim)) 164 else if (maxwell_field == OPTION__USERDEFINEDINITIALMAXWELLSTATES__MAGNETIC_FIELD) then 165 call parse_block_string( blk, il - 1, 3, st%user_def_b_field(idim)) 166 call messages_write(" B-field in dimension "//trim(cdim)//" : "//trim(st%user_def_b_field(idim)), fmt='(a,i1,2a)') 167 call conv_to_C_string(st%user_def_b_field(idim)) 168 end if 169 ! fill Maxwell states with user-defined formulas 170 do ip = 1, mesh%np_part 171 xx = mesh%x(ip, :) 172 rr = sqrt(sum(xx(:)**2)) 173 ! parse user-defined expressions 174 if (maxwell_field == OPTION__USERDEFINEDINITIALMAXWELLSTATES__ELECTRIC_FIELD) then 175 call parse_expression(e_value, dummy, st%dim, xx, rr, M_ZERO, & 176 st%user_def_e_field(idim)) 177 b_value = M_ZERO 178 else if (maxwell_field == OPTION__USERDEFINEDINITIALMAXWELLSTATES__MAGNETIC_FIELD) then 179 call parse_expression(b_value, dummy, st%dim, xx, rr, M_ZERO, & 180 st%user_def_b_field(idim)) 181 e_value = M_ZERO 182 end if 183 e_value = units_to_atomic(units_inp%energy/units_inp%length, e_value) 184 b_value = units_to_atomic(unit_one/(units_inp%length**2), b_value) 185 ! fill state 186 call build_rs_element(e_value, b_value, st%rs_sign, rs_state_add(ip,idim), & 187 st%ep(ip), st%mu(ip)) 188 end do 189 190 case(STATE_FROM_FILE) 191 ! The input format can be coded in column four now. As it is 192 ! not used now, we just say "file". 193 ! Read the filename. 194 call parse_block_integer(blk, il - 1, 2, maxwell_field ) 195 SAFE_ALLOCATE(e_field(1:mesh%np)) 196 SAFE_ALLOCATE(b_field(1:mesh%np)) 197 e_field = M_ZERO 198 b_field = M_ZERO 199 if (maxwell_field == OPTION__USERDEFINEDINITIALMAXWELLSTATES__ELECTRIC_FIELD) then 200 call parse_block_string(blk, il - 1, 3, filename_e_field) 201 call messages_write(" E-field in dimension "//trim(cdim)//" : "//trim(filename_e_field), fmt='(a,i1,2a)') 202 call dio_function_input(filename_e_field, namespace, mesh, e_field(:), ierr) 203 if (ierr > 0) then 204 message(1) = 'Could not read the file!' 205 write(message(2),'(a,i1)') 'Error code: ', ierr 206 call messages_fatal(2, namespace=namespace) 207 end if 208 e_field = units_to_atomic(units_inp%energy/units_inp%length, e_field) 209 else if (maxwell_field == OPTION__USERDEFINEDINITIALMAXWELLSTATES__MAGNETIC_FIELD) then 210 call parse_block_string(blk, il - 1, 3, filename_b_field) 211 call messages_write(" B-field in dimension "//trim(cdim)//" : "//trim(filename_b_field), fmt='(a,i1,2a)') 212 call dio_function_input(filename_b_field, namespace, mesh, b_field(:), ierr) 213 if (ierr > 0) then 214 message(1) = 'Could not read the file!' 215 write(message(2),'(a,i1)') 'Error code: ', ierr 216 call messages_fatal(2, namespace=namespace) 217 end if 218 b_field = units_to_atomic(unit_one/units_inp%length**2, b_field) 219 end if 220 ! fill state 221 call build_rs_vector(e_field(:), b_field(:), st%rs_sign, rs_state_add(:,idim), & 222 st%ep(ip), st%mu(ip)) 223 224 SAFE_DEALLOCATE_A(e_field) 225 SAFE_DEALLOCATE_A(b_field) 226 227 case default 228 message(1) = 'Wrong entry in UserDefinedMaxwellStates, column 2.' 229 message(2) = 'You may state "formula" or "file" here.' 230 call messages_fatal(2, namespace=namespace) 231 end select 232 233 rs_state(:,idim) = rs_state(:,idim) + rs_state_add(:,idim) 234 235 end do 236 237 do ip = 1, mesh%np 238 user_def_rs_state(ip,:) = rs_state(ip,:) 239 end do 240 241 SAFE_DEALLOCATE_A(rs_state) 242 SAFE_DEALLOCATE_A(rs_state_add) 243 call parse_block_end(blk) 244 !call messages_print_stress(stdout, namespace=namespace) 245 246 else 247 message(1) = "'UserDefineInitialdStates' has to be specified as block." 248 call messages_fatal(1, namespace=namespace) 249 end if 250 251 POP_SUB(states_mxll_read_user_def) 252 end subroutine states_mxll_read_user_def 253 254 255 !---------------------------------------------------------- 256 subroutine states_mxll_dump(restart, st, gr, zff, zff_dim, ierr, iter, st_start_writing, verbose) 257 type(restart_t), intent(in) :: restart 258 type(states_mxll_t), intent(in) :: st 259 type(grid_t), intent(in) :: gr 260 CMPLX, intent(in) :: zff(:,:) 261 integer, intent(in) :: zff_dim 262 integer, intent(out) :: ierr 263 integer, optional, intent(in) :: iter 264 integer, optional, intent(in) :: st_start_writing 265 logical, optional, intent(in) :: verbose 266 267 integer :: iunit_wfns, iunit_states 268 integer :: err, err2(2), ist, idim, itot 269 integer :: root(1:P_STRATEGY_MAX) 270 character(len=MAX_PATH_LEN) :: filename 271 character(len=300) :: lines(3) 272 logical :: should_write, verbose_ 273 274 PUSH_SUB(states_mxll_dump) 275 276 verbose_ = optional_default(verbose, .true.) 277 278 ierr = 0 279 280 if(restart_skip(restart)) then 281 POP_SUB(states_mxll_dump) 282 return 283 end if 284 285 if(verbose_) then 286 message(1) = "Info: Writing Maxwell states." 287 call print_date(trim(message(1))//' ') 288 end if 289 290 call profiling_in(prof_write, "MAXWELL_RESTART_WRITE") 291 292 call restart_block_signals() 293 294 iunit_states = restart_open(restart, 'maxwell_states') 295 write(lines(2), '(a20,1i10)') 'dim= ', zff_dim 296 call restart_write(restart, iunit_states, lines, 3, err) 297 if (err /= 0) ierr = ierr + 1 298 call restart_close(restart, iunit_states) 299 300 iunit_wfns = restart_open(restart, 'wfns') 301 lines(1) = '# #dim filename' 302 lines(2) = '%RS States' 303 call restart_write(restart, iunit_wfns, lines, 2, err) 304 if (err /= 0) ierr = ierr + 2 305 306 307 itot = 1 308 root = 0 309 err2 = 0 310 ist = 1 311 312 do idim = 1, zff_dim 313 itot = itot + 1 314 315 root(P_STRATEGY_DOMAINS) = mod(itot - 1, gr%mesh%mpi_grp%size) 316 write(filename,'(i10.10)') itot 317 318 write(lines(1), '(i8,3a)') idim, ' | "', trim(filename), '"' 319 call restart_write(restart, iunit_wfns, lines, 1, err) 320 if (err /= 0) err2(1) = err2(1) + 1 321 322 should_write = st%st_start <= ist .and. ist <= st%st_end 323 if (should_write .and. present(st_start_writing)) then 324 if (ist < st_start_writing) should_write = .false. 325 end if 326 327 if (should_write) then 328 call zrestart_write_mesh_function(restart, filename, gr%mesh, zff(:,idim), err, root = root) 329 if (err /= 0) err2(2) = err2(2) + 1 330 end if 331 end do ! zff_dim 332 333 if (err2(1) /= 0) ierr = ierr + 8 334 if (err2(2) /= 0) ierr = ierr + 16 335 336 lines(1) = '%' 337 call restart_write(restart, iunit_wfns, lines, 1, err) 338 if (err /= 0) ierr = ierr + 64 339 if (present(iter)) then 340 write(lines(1),'(a,i7)') 'Iter = ', iter 341 call restart_write(restart, iunit_wfns, lines, 1, err) 342 if (err /= 0) ierr = ierr + 128 343 end if 344 345 call restart_close(restart, iunit_wfns) 346 347 if(verbose_) then 348 message(1) = "Info: Finished writing Maxwell states." 349 call print_date(trim(message(1))//' ') 350 end if 351 352 call restart_unblock_signals() 353 354 call profiling_out(prof_write) 355 POP_SUB(states_mxll_dump) 356 return 357 358 end subroutine states_mxll_dump 359 360 !---------------------------------------------------------- 361 subroutine states_mxll_load(restart, st, gr, namespace, zff, zff_dim, ierr, iter, lowest_missing, label, verbose) 362 type(restart_t), intent(in) :: restart 363 type(states_mxll_t), intent(inout) :: st 364 type(grid_t), intent(in) :: gr 365 type(namespace_t), intent(in) :: namespace 366 CMPLX, intent(inout) :: zff(:,:) 367 integer, intent(in) :: zff_dim 368 integer, intent(out) :: ierr 369 integer, optional, intent(out) :: iter 370 integer, optional, intent(out) :: lowest_missing(:) !< all states below this one were read successfully 371 character(len=*), optional, intent(in) :: label 372 logical, optional, intent(in) :: verbose 373 374 integer :: states_file, wfns_file, err, ist, idim, dim, mx_st_start, mx_st_end 375 integer :: idone, iread, ntodo 376 character(len=12) :: filename 377 character(len=1) :: char 378 logical, allocatable :: filled(:, :) 379 character(len=256) :: lines(3), label_ 380 character(len=50) :: str 381 382 logical :: verbose_ 383 character(len=256), allocatable :: restart_file(:, :) 384 logical, allocatable :: restart_file_present(:, :) 385 386 PUSH_SUB(states_mxll_load) 387 388 ierr = 0 389 dim = zff_dim 390 391 ! make sure these intent(out)`s are initialized no matter what 392 if (present(lowest_missing)) lowest_missing = 1 393 if (present(iter)) iter = 0 394 395 if (restart_skip(restart)) then 396 ierr = -1 397 POP_SUB(states_load) 398 return 399 end if 400 401 call profiling_in(prof_read, "RESTART_READ") 402 403 verbose_ = optional_default(verbose, .true.) 404 405 if (present(label)) then 406 label_ = trim(label) 407 end if 408 409 message(1) = 'Info: Reading Maxwell states' 410 if (len(trim(label_)) > 0) then 411 message(1) = trim(message(1)) // trim(label_) 412 end if 413 message(1) = trim(message(1)) // "." 414 if(verbose_) call print_date(trim(message(1))//' ') 415 416 states_file = restart_open(restart, 'maxwell_states') 417 call restart_read(restart, states_file, lines, 3, err) 418 if (err /= 0) then 419 ierr = ierr - 2 420 else 421 read(lines(2), *) idim 422 end if 423 call restart_close(restart, states_file) 424 425 ! open files to read 426 wfns_file = restart_open(restart, 'wfns') 427 call restart_read(restart, wfns_file, lines, 2, err) 428 if (err /= 0) then 429 ierr = ierr - 2**5 430 end if 431 432 ! If any error occured up to this point then it is not worth continuing, 433 ! as there something fundamentally wrong with the restart files 434 if (ierr /= 0) then 435 call restart_close(restart, wfns_file) 436 call profiling_out(prof_read) 437 POP_SUB(states_mxll_load) 438 return 439 end if 440 441 SAFE_ALLOCATE(restart_file(1:zff_dim, st%st_start:st%st_end)) 442 SAFE_ALLOCATE(restart_file_present(1:zff_dim,st%st_start:st%st_end)) 443 restart_file_present = .false. 444 445 ! Next we read the list of states from the files. 446 ! Errors in reading the information of a specific state from the files are ignored 447 ! at this point, because later we will skip reading the wavefunction of that state. 448 do 449 call restart_read(restart, wfns_file, lines, 1, err) 450 if (err == 0) then 451 read(lines(1), '(a)') char 452 if (char == '%') then 453 !We reached the end of the file 454 exit 455 else 456 read(lines(1), *) idim, char, filename 457 end if 458 end if 459 460 if (ist >= st%st_start .and. ist <= st%st_end) then 461 restart_file(idim, ist) = trim(filename) 462 restart_file_present(idim, ist) = .true. 463 end if 464 465 end do 466 if (present(iter)) then 467 call restart_read(restart, wfns_file, lines, 1, err) 468 if (err /= 0) then 469 ierr = ierr - 2**8 470 else 471 read(lines(1), *) filename, filename, iter 472 end if 473 end if 474 475 call restart_close(restart, wfns_file) 476 477 ! Now we read the wavefunctions. At this point we need to have all the information from the 478 ! states, and wfns files in order to avoid serialisation of reads, as restart_read 479 ! works as a barrier. 480 481 mx_st_start=st%st_start 482 mx_st_end=st%st_end 483 SAFE_ALLOCATE(filled(1:zff_dim,mx_st_start:mx_st_end)) 484 filled = .false. 485 486 if (present(lowest_missing)) lowest_missing = st%nst + 1 487 488 iread = 0 489 if (mpi_grp_is_root(mpi_world) .and. verbose_) then 490 idone = 0 491 ntodo = st%lnst*zff_dim 492 call loct_progress_bar(-1, ntodo) 493 end if 494 495 ist = 1 496 do idim = 1, zff_dim 497 if (.not. restart_file_present(idim, ist)) then 498 if (present(lowest_missing)) & 499 lowest_missing(idim) = min(lowest_missing(idim), ist) 500 cycle 501 endif 502 503 call zrestart_read_mesh_function(restart, restart_file(idim, ist), gr%mesh, & 504 zff(:,idim), err) 505 506 507 if (err == 0) then 508 filled(idim, ist) = .true. 509 iread = iread + 1 510 else if (present(lowest_missing)) then 511 lowest_missing(idim) = min(lowest_missing(idim), ist) 512 end if 513 514 if (mpi_grp_is_root(mpi_world) .and. verbose_) then 515 idone = idone + 1 516 call loct_progress_bar(idone, ntodo) 517 end if 518 519 end do 520 521 SAFE_DEALLOCATE_A(restart_file) 522 SAFE_DEALLOCATE_A(restart_file_present) 523 SAFE_DEALLOCATE_A(filled) 524 525 if(mpi_grp_is_root(mpi_world) .and. verbose_) then 526 call messages_new_line() 527 end if 528 529 if (ierr == 0 .and. iread /= st%nst * zff_dim) then 530 if(iread > 0) then 531 ierr = iread 532 else 533 ierr = -1 534 endif 535 ! otherwise ierr = 0 would mean either all was read correctly, or nothing at all was read! 536 537 write(str, '(a,i5)') 'Reading Maxwell states.' 538 call messages_print_stress(stdout, trim(str), namespace=namespace) 539 write(message(1),'(a,i6,a,i6,a)') 'Only ', iread,' files out of ', & 540 st%nst * zff_dim, ' could be read.' 541 call messages_info(1) 542 call messages_print_stress(stdout, namespace=namespace) 543 end if 544 545 message(1) = 'Info: Maxwell states reading done.' 546 if(verbose_) call print_date(trim(message(1))//' ') 547 548 call profiling_out(prof_read) 549 POP_SUB(states_mxll_load) 550 551 end subroutine states_mxll_load 552 553end module states_mxll_restart_oct_m 554 555 556!! Local Variables: 557!! mode: f90 558!! coding: utf-8 559!! End: 560