Octopus
states_elec_restart.F90
Go to the documentation of this file.
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
22
24#ifdef HAVE_ADIOS2
25 use adios2
26#endif
27 use batch_oct_m
29 use debug_oct_m
31 use global_oct_m
33 use, intrinsic :: iso_fortran_env
37 use loct_oct_m
38 use math_oct_m
39 use mesh_oct_m
43 use mpi_oct_m
47 use parser_oct_m
51 use smear_oct_m
52 use space_oct_m
57 use string_oct_m
58 use types_oct_m
60
61 implicit none
62
63
64 private
65 public :: &
77
78contains
79
80 ! ---------------------------------------------------------
81 subroutine states_elec_look_and_load(restart, namespace, space, st, mesh, kpoints, is_complex, packed)
82 type(restart_t), intent(in) :: restart
83 type(namespace_t), intent(in) :: namespace
84 class(space_t), intent(in) :: space
85 type(states_elec_t), target, intent(inout) :: st
86 class(mesh_t), intent(in) :: mesh
87 type(kpoints_t), intent(in) :: kpoints
88 logical, optional, intent(in) :: is_complex
89 logical, optional, intent(in) :: packed
90
91 integer :: nkpt, dim, nst, ierr
92 real(real64), allocatable :: new_occ(:,:)
93
95
96 !check how many wfs we have
97 call states_elec_look(restart, nkpt, dim, nst, ierr)
98 if (ierr /= 0) then
99 message(1) = "Unable to read states information."
100 call messages_fatal(1, namespace=namespace)
101 end if
102
103 if (st%parallel_in_states) then
104 message(1) = "Internal error: cannot use states_elec_look_and_load when parallel in states."
105 call messages_fatal(1, namespace=namespace)
106 end if
107
108 ! Resize st%occ, retaining information there
109 allocate(new_occ(1:nst, 1:st%nik))
110 new_occ(:,:) = m_zero
111 new_occ(1:min(nst, st%nst),:) = st%occ(1:min(nst, st%nst),:)
112 safe_deallocate_a(st%occ)
113 call move_alloc(new_occ, st%occ)
114
115 ! FIXME: This wrong, one cannot just change the number of states
116 ! without updating the internal structures, in the case of parallel in states.
117 ! I think it is right now without state parallelism.
118 st%nst = nst
119 st%st_start = 1
120 st%st_end = nst
121 st%lnst = nst
122
123 safe_deallocate_a(st%node)
124 safe_allocate(st%node(1:st%nst))
125 st%node(:) = 0
126
127 safe_deallocate_a(st%eigenval)
128 safe_allocate(st%eigenval(1:st%nst, 1:st%nik))
129 st%eigenval = huge(st%eigenval)
130
131 if (present(is_complex)) then
132 if (is_complex) then
133 call states_elec_allocate_wfns(st, mesh, type_cmplx, packed=optional_default(packed, .false.))
134 else
135 call states_elec_allocate_wfns(st, mesh, type_float, packed=optional_default(packed, .false.))
136 end if
137 else
138 ! allow states_elec_allocate_wfns to decide for itself whether complex or real needed
139 call states_elec_allocate_wfns(st, mesh, packed=optional_default(packed, .false.))
140 end if
141
142 if (st%d%ispin == spinors) then
143 safe_allocate(st%spin(1:3, 1:st%nst, 1:st%nik))
144 st%spin = m_zero
145 end if
146
147 ! load wavefunctions
148 call states_elec_load(restart, namespace, space, st, mesh, kpoints, ierr)
149 if (ierr /= 0) then
150 message(1) = "Unable to read wavefunctions."
151 call messages_fatal(1, namespace=namespace)
152 end if
153
155 end subroutine states_elec_look_and_load
156
157
158 ! ---------------------------------------------------------
159 subroutine states_elec_dump(restart, space, st, mesh, kpoints, ierr, iter, lr, verbose)
160 type(restart_t), intent(in) :: restart
161 class(space_t), intent(in) :: space
162 type(states_elec_t), target, intent(in) :: st
163 class(mesh_t), intent(in) :: mesh
164 type(kpoints_t), intent(in) :: kpoints
165 integer, intent(out) :: ierr
166 integer, optional, intent(in) :: iter
167 type(lr_t), optional, intent(in) :: lr
168 logical, optional, intent(in) :: verbose
169
170 integer :: iunit_wfns, iunit_occs, iunit_states
171 integer :: err, err2(2), ik, idir, ist, idim, itot
172 integer :: root(1:P_STRATEGY_MAX)
173 character(len=MAX_PATH_LEN) :: filename
174 character(len=300) :: lines(3)
175 logical :: lr_wfns_are_associated, should_write, verbose_
176 real(real64) :: kpoint(space%dim)
177 real(real64), allocatable :: dpsi(:), rff_global(:)
178 complex(real64), allocatable :: zpsi(:), zff_global(:)
179 integer :: restart_file_format
180
181 push_sub(states_elec_dump)
182
183 verbose_ = optional_default(verbose, .true.)
184
185 ierr = 0
186
187 if (restart%skip()) then
188 pop_sub(states_elec_dump)
189 return
190 end if
191
192 if (verbose_) then
193 message(1) = "Info: Writing states."
194 call print_date(trim(message(1))//' ')
195 end if
196
197 call profiling_in("RESTART_WRITE")
198
199 restart_file_format = restart%file_format_states
200 if (present(lr)) then
201 lr_wfns_are_associated = (allocated(lr%ddl_psi) .and. states_are_real(st)) .or. &
202 (allocated(lr%zdl_psi) .and. states_are_complex(st))
203 assert(lr_wfns_are_associated)
204 ! always use obf format for linear response
205 restart_file_format = option__restartfileformatstates__obf
206 end if
207
209
210 iunit_states = restart%open('states')
211 write(lines(1), '(a20,1i10)') 'nst= ', st%nst
212 write(lines(2), '(a20,1i10)') 'dim= ', st%d%dim
213 write(lines(3), '(a20,1i10)') 'nik= ', st%nik
214 call restart%write(iunit_states, lines, 3, err)
215 if (err /= 0) ierr = ierr + 1
216 call restart%close(iunit_states)
217
218
219 iunit_wfns = restart%open('wfns')
220 lines(1) = '# #k-point #st #dim filename'
221 if (states_are_real(st)) then
222 lines(2) = '%Real_Wavefunctions'
223 else
224 lines(2) = '%Complex_Wavefunctions'
225 end if
226 call restart%write(iunit_wfns, lines, 2, err)
227 if (err /= 0) ierr = ierr + 2
228
229
230 iunit_occs = restart%open('occs')
231 lines(1) = '# occupations | eigenvalue[a.u.] | Im(eigenvalue) [a.u.] | k-points | k-weights | filename | ik | ist | idim'
232 lines(2) = '%Occupations_Eigenvalues_K-Points'
233 call restart%write(iunit_occs, lines, 2, err)
234 if (err /= 0) ierr = ierr + 4
235
236
237 if (states_are_real(st)) then
238 safe_allocate(dpsi(1:mesh%np))
239 safe_allocate(rff_global(1:mesh%np_global))
240 else
241 safe_allocate(zpsi(1:mesh%np))
242 safe_allocate(zff_global(1:mesh%np_global))
243 end if
244
245 itot = 1
246 root = -1
247 err2 = 0
248 do ik = 1, st%nik
249 kpoint(1:space%dim) = &
250 kpoints%get_point(st%d%get_kpoint_index(ik), absolute_coordinates = .true.)
251
252 do ist = 1, st%nst
253 do idim = 1, st%d%dim
254 root(p_strategy_domains) = mod(itot - 1, mesh%mpi_grp%size)
255 write(filename,'(i10.10)') itot
256
257 write(lines(1), '(i8,a,i8,a,i8,3a)') ik, ' | ', ist, ' | ', idim, ' | "', trim(filename), '"'
258 call restart%write(iunit_wfns, lines, 1, err)
259 if (err /= 0) err2(1) = err2(1) + 1
260
261 write(lines(1), '(e23.16,a,e23.16)') st%occ(ist,ik), ' | ', st%eigenval(ist, ik)
262 write(lines(1), '(a,a,e23.16)') trim(lines(1)), ' | ', m_zero
263 do idir = 1, space%dim
264 write(lines(1), '(a,a,e23.16)') trim(lines(1)), ' | ', kpoint(idir)
265 end do
266 write(lines(1), '(a,a,e23.16,a,i10.10,3(a,i8))') trim(lines(1)), &
267 ' | ', st%kweights(ik), ' | ', itot, ' | ', ik, ' | ', ist, ' | ', idim
268 call restart%write(iunit_occs, lines, 1, err)
269 if (err /= 0) err2(1) = err2(1) + 1
270
271 should_write = st%st_start <= ist .and. ist <= st%st_end
272
273 if (restart_file_format == option__restartfileformatstates__obf) then
274 if (should_write) then
275 if (st%d%kpt%start <= ik .and. ik <= st%d%kpt%end) then
276 if (states_are_real(st)) then
277 if (.not. present(lr)) then
278 call states_elec_get_state(st, mesh, idim, ist, ik, dpsi)
279 call restart%write_mesh_function(filename, mesh, dpsi, err, root = root)
280 else
281 call restart%write_mesh_function(filename, mesh, &
282 lr%ddl_psi(:, idim, ist, ik), err, root = root)
283 end if
284 else
285 if (.not. present(lr)) then
286 call states_elec_get_state(st, mesh, idim, ist, ik, zpsi)
287 call restart%write_mesh_function(filename, mesh, zpsi, err, root = root)
288 else
289 call restart%write_mesh_function(filename, mesh, &
290 lr%zdl_psi(:, idim, ist, ik), err, root = root)
291 end if
292 end if
293 else
294 err = 0
295 end if
296 end if
297 if (err /= 0) err2(2) = err2(2) + 1
298 end if
299
300 itot = itot + 1
301 end do ! st%d%dim
302 end do ! st%nst
303 end do ! st%nik
304 if (err2(1) /= 0) ierr = ierr + 8
305 if (err2(2) /= 0) ierr = ierr + 16
306
307 safe_deallocate_a(dpsi)
308 safe_deallocate_a(zpsi)
309 safe_deallocate_a(rff_global)
310 safe_deallocate_a(zff_global)
311
312 if (restart_file_format == option__restartfileformatstates__adios2) then
313 if (present(lr)) then
314 ! this should never occur because we force obf format for linear response
315 call messages_not_implemented("ADIOS2 restart file format with linear response.")
316 end if
317 call states_elec_dump_adios2(restart, st, mesh, ierr)
318 end if
319
320 lines(1) = '%'
321 call restart%write(iunit_occs, lines, 1, err)
322 if (err /= 0) ierr = ierr + 32
323 call restart%write(iunit_wfns, lines, 1, err)
324 if (err /= 0) ierr = ierr + 64
325 if (present(iter)) then
326 write(lines(1),'(a,i7)') 'Iter = ', iter
327 call restart%write(iunit_wfns, lines, 1, err)
328 if (err /= 0) ierr = ierr + 128
329 end if
330
331 call restart%close(iunit_wfns)
332 call restart%close(iunit_occs)
333
334 if (verbose_) then
335 message(1) = "Info: Finished writing states."
336 call print_date(trim(message(1))//' ')
337 end if
338
340
341 call profiling_out("RESTART_WRITE")
342 pop_sub(states_elec_dump)
343 end subroutine states_elec_dump
344
345 ! ---------------------------------------------------------
346 subroutine states_elec_dump_adios2(restart, st, mesh, ierr)
347 type(restart_t), intent(in) :: restart
348 type(states_elec_t), target, intent(in) :: st
349 class(mesh_t), intent(in) :: mesh
350 integer, intent(out) :: ierr
351
352#ifdef HAVE_ADIOS2
353 type(adios2_adios) :: adios
354 type(adios2_io) :: io
355 type(adios2_variable) :: var, var_indices
356 type(adios2_attribute) :: attribute
357 type(adios2_engine) :: engine
358 integer :: adios2_type, ik, ib, ip, adios2_mode
359 integer(int64), allocatable :: global_indices(:)
360 class(wfs_elec_t), pointer :: psib
361#endif
363
364#ifdef HAVE_ADIOS2
365#ifdef HAVE_MPI
366 call adios2_init(adios, restart%mpi_grp%comm%MPI_VAL, ierr)
367#else
368 call adios2_init(adios, ierr)
369#endif
370 call check_error(.not. adios%valid .or. ierr /= adios2_error_none, "Problem initializing ADIOS2.")
371 call adios2_declare_io(io, adios, "writer", ierr)
372 call check_error(.not. io%valid .or. ierr /= adios2_error_none, "Problem initializing ADIOS2.")
373 ! do not compute min and max
374 call adios2_set_parameter(io, "StatsLevel", "0", ierr)
375 call check_error(ierr /= adios2_error_none, "Problem setting ADIOS2 parameter.")
376 ! open restart file
377 call adios2_open(engine, io, trim(restart%dir())//"/"//"wavefunctions.bp", adios2_mode_write, ierr)
378 call check_error(.not. engine%valid .or. ierr /= adios2_error_none, "Problem opening ADIOS2 restart file.")
379 if (states_are_real(st)) then
380 adios2_type = adios2_type_dp
381 else
382 adios2_type = adios2_type_complex_dp
383 end if
384 ! define one variable for all states
385 call adios2_define_variable(var, io, "wavefunctions", adios2_type, 3, &
386 [int(st%nst * st%d%dim, int64), int(mesh%np_global, int64), int(st%nik, int64)], & ! global size of wavefunctions
387 [0_int64, 0_int64, 0_int64], & ! local selection of the variable, overwritten by set_selection call below
388 [1_int64, 1_int64, 1_int64], &
389 adios2_variable_dims, ierr)
390 call check_error(.not. var%valid .or. ierr /= adios2_error_none, "Problem creating ADIOS2 variable (wavefunctions).")
391 ! save ParDomains in an attribute
392 call adios2_define_attribute(attribute, io, "ParDomains", mesh%mpi_grp%size, ierr)
393 call check_error(.not. attribute%valid .or. ierr /= adios2_error_none, "Problem creating ADIOS2 attribute (ParDomains).")
394 ! define one variable for all states
395 call adios2_define_variable(var_indices, io, "global_indices", adios2_type_integer8, 1, &
396 [int(mesh%np_global, int64)], & ! global size of array
397 [int(mesh%pv%xlocal-1, int64)], & ! local offset
398 [int(mesh%np, int64)], & ! local size
399 adios2_variable_dims, ierr)
400 call check_error(.not. var_indices%valid .or. ierr /= adios2_error_none, "Problem creating ADIOS2 variable (global_indices).")
401 ! save grid indices
402 safe_allocate(global_indices(mesh%np))
403 do ip = 1, mesh%np
404 global_indices(ip) = par_vec_local2global(mesh%pv, ip)
405 end do
406
407 call adios2_begin_step(engine, ierr)
408 call adios2_put(engine, var_indices, global_indices, ierr)
409
410 do ik = st%d%kpt%start, st%d%kpt%end
411 do ib = st%group%block_start, st%group%block_end
412 ! make sure we have a batch in BATCH_PACKED status
413 select case (st%group%psib(ib, ik)%status())
414 case (batch_not_packed)
415 safe_allocate(psib)
416 call st%group%psib(ib, ik)%copy_to(psib)
417 ! transform the unpacked batch to a packed batch on the CPU
418 call psib%do_pack(batch_packed)
419 ! we need to use sync mode in this case because we use a temporary batch
420 adios2_mode = adios2_mode_sync
421 case (batch_packed)
422 psib => st%group%psib(ib, ik)
423 adios2_mode = adios2_mode_deferred
425 safe_allocate(psib)
426 call st%group%psib(ib, ik)%copy_to(psib, copy_data=.true.)
427 ! copy to CPU
428 call psib%do_unpack(force=.true.)
429 ! make sure the batch is in packed state on the CPU
430 call psib%do_pack(batch_packed)
431 ! we need to use sync mode in this case because we use a temporary batch
432 adios2_mode = adios2_mode_sync
433 end select
434 assert(psib%status() == batch_packed)
435
436 ! select part of variable to write, correspond to the memory of the current batch
437 ! first argument gives the start in the global array (-1 needed for C indexing in ADIOS2)
438 ! second argument gives the size of the local batch
439 call adios2_set_selection(var, 3, &
440 [int((states_elec_block_min(st, ib)-1)*st%d%dim, int64), int(mesh%pv%xlocal-1, int64), int(ik-1, int64)], &
441 [int(st%group%psib(ib, ik)%nst_linear, int64), int(mesh%np, int64), 1_int64], ierr)
442 ! set memory layout of one batch, take into account padding and ghost/boundary points
443 if (states_are_real(st)) then
444 call adios2_set_memory_selection(var, 3, [0_int64, 0_int64, 0_int64], &
445 [int(size(psib%dff_pack, 1), int64), &
446 int(size(psib%dff_pack, 2), int64), 1_int64], ierr)
447 call adios2_put(engine, var, psib%dff_pack, adios2_mode, ierr)
448 else
449 call adios2_set_memory_selection(var, 3, [0_int64, 0_int64, 0_int64], &
450 [int(size(psib%zff_pack, 1), int64), &
451 int(size(psib%zff_pack, 2), int64), 1_int64], ierr)
452 call adios2_put(engine, var, psib%zff_pack, adios2_mode, ierr)
453 end if
454 select case (st%group%psib(ib, ik)%status())
456 call psib%end()
457 safe_deallocate_p(psib)
458 case (batch_packed)
459 nullify(psib)
460 end select
461 end do
462 end do
463 call adios2_end_step(engine, ierr)
464
465 ! deallocate only after end_step because only now everything has been written
466 safe_deallocate_a(global_indices)
467
468 ! close restart file
469 call adios2_close(engine, ierr)
470 call check_error(ierr /= adios2_error_none, "Problem closing ADIOS2 engine.")
471 call adios2_finalize(adios, ierr)
472 call check_error(ierr /= adios2_error_none, "Problem finalizing ADIOS2.")
473#else
474 ! avoid empty routine
475 ierr = 0
476#endif
477
479 end subroutine states_elec_dump_adios2
480
481
482 ! ---------------------------------------------------------
487 subroutine states_elec_load(restart, namespace, space, st, mesh, kpoints, ierr, iter, lr, lowest_missing, label, verbose, skip)
488 type(restart_t), intent(in) :: restart
489 type(namespace_t), intent(in) :: namespace
490 class(space_t), intent(in) :: space
491 type(states_elec_t), target, intent(inout) :: st
492 class(mesh_t), intent(in) :: mesh
493 type(kpoints_t), intent(in) :: kpoints
494 integer, intent(out) :: ierr
495 integer, optional, intent(out) :: iter
496 type(lr_t), optional, intent(inout) :: lr
497 integer, optional, intent(out) :: lowest_missing(:, :)
498 character(len=*), optional, intent(in) :: label
499 logical, optional, intent(in) :: verbose
500 logical, optional, intent(in) :: skip(:)
501
502 integer :: states_elec_file, wfns_file, occ_file, err, ik, ist, idir, idim
503 integer :: idone, iread, ntodo
504 character(len=12) :: filename
505 character(len=1) :: char
506 logical, allocatable :: filled(:, :, :)
507 character(len=256) :: lines(3), label_
508 character(len=50) :: str
509
510 real(real64) :: my_occ, imev, my_kweight
511 logical :: read_occ, lr_allocated, verbose_
512 logical :: integral_occs
513 real(real64), allocatable :: dpsi(:)
514 complex(real64), allocatable :: zpsi(:), zpsil(:)
515 character(len=256), allocatable :: restart_file(:, :, :)
516 logical, allocatable :: restart_file_present(:, :, :)
517 real(real64) :: kpoint(space%dim), read_kpoint(space%dim)
518
519 integer :: iread_tmp
520 integer, allocatable :: lowest_missing_tmp(:, :)
521 integer :: restart_file_format
522
523 push_sub(states_elec_load)
524
525 ierr = 0
526
527 ! make sure these intent(out)`s are initialized no matter what
528 if (present(lowest_missing)) lowest_missing = 1
529 if (present(iter)) iter = 0
530
531 if (present(skip)) then
532 assert(ubound(skip, dim = 1) == st%nst)
533 end if
534
535 if (restart%skip()) then
536 ierr = -1
537 pop_sub(states_elec_load)
538 return
539 end if
540
541 call profiling_in("RESTART_READ")
542
543 verbose_ = optional_default(verbose, .true.)
544
545 if (present(label)) then
546 label_ = trim(label)
547 else
548 if (present(lr)) then
549 label_ = " for linear response"
550 else
551 label_ = ""
552 end if
553 end if
554
555 message(1) = 'Info: Reading states'
556 if (len(trim(label_)) > 0) then
557 message(1) = trim(message(1)) // trim(label_)
558 end if
559 message(1) = trim(message(1)) // "."
560 if (verbose_) call print_date(trim(message(1))//' ')
561
562 if (.not. present(lr)) then
563 st%fromScratch = .false. ! obviously, we are using restart info
564 end if
565
566 ! If one restarts a GS calculation changing the %Occupations block, one
567 ! cannot read the occupations, otherwise these overwrite the ones from
568 ! the input file. restart_fixed_occ makes that we do use the ones in the file.
569 integral_occs = .true. ! only used if restart_fixed_occ
570 if (st%restart_fixed_occ) then
571 read_occ = .true.
572 st%fixed_occ = .true.
573 else
574 read_occ = .not. st%fixed_occ
575 end if
576
577 if (.not. present(lr)) then
578 st%eigenval(:, :) = m_zero
579 ! to be filled in from reading afterward
580 end if
581
582 if (.not. present(lr) .and. read_occ) then
583 st%occ(:, :) = m_zero
584 ! to be filled in from reading afterward
585 end if
586
587 restart_file_format = restart%file_format_states
588 ! sanity check
589 if (present(lr)) then
590 lr_allocated = (allocated(lr%ddl_psi) .and. states_are_real(st)) .or. &
591 (allocated(lr%zdl_psi) .and. states_are_complex(st))
592 assert(lr_allocated)
593 ! always use obf format for linear response
594 restart_file_format = option__restartfileformatstates__obf
595 end if
596
597 states_elec_file = restart%open('states')
598 ! sanity check on spin/k-points. Example file 'states':
599 ! nst= 2
600 ! dim= 1
601 ! nik= 2
602 call restart%read(states_elec_file, lines, 3, err)
603 if (err /= 0) then
604 ierr = ierr - 2
605 else
606 read(lines(2), *) str, idim
607 read(lines(3), *) str, ik
608 if (idim == 2 .and. st%d%dim == 1) then
609 write(message(1),'(a)') 'Incompatible restart information: saved calculation is spinors, this one is not.'
610 call messages_warning(1, namespace=namespace)
611 ierr = ierr - 2**2
612 end if
613 if (idim == 1 .and. st%d%dim == 2) then
614 write(message(1),'(a)') 'Incompatible restart information: this calculation is spinors, saved one is not.'
615 call messages_warning(1, namespace=namespace)
616 ierr = ierr - 2**3
617 end if
618 if (ik < st%nik) then
619 write(message(1),'(a)') 'Incompatible restart information: not enough k-points.'
620 write(message(2),'(2(a,i6))') 'Expected ', st%nik, ' > Read ', ik
621 call messages_warning(2, namespace=namespace)
622 end if
623 ! We will check that they are the right k-points later, so we do not need to put a specific error here.
624 end if
625 call restart%close(states_elec_file)
626
627
628 ! open files to read
629 wfns_file = restart%open('wfns')
630 occ_file = restart%open('occs')
631 call restart%read(wfns_file, lines, 2, err)
632 if (err /= 0) then
633 ierr = ierr - 2**5
634 else if (states_are_real(st)) then
635 read(lines(2), '(a)') str
636 if (str(2:8) == 'Complex') then
637 message(1) = "Cannot read real states from complex wavefunctions."
638 call messages_warning(1, namespace=namespace)
639 ierr = ierr - 2**6
640 else if (str(2:5) /= 'Real') then
641 message(1) = "Restart file 'wfns' does not specify real/complex; cannot check compatibility."
642 call messages_warning(1, namespace=namespace)
643 end if
644 end if
645 ! complex can be restarted from real, so there is no problem.
646
647 ! Skip two lines.
648 call restart%read(occ_file, lines, 2, err)
649 if (err /= 0) ierr = ierr - 2**7
650
651 ! If any error occured up to this point then it is not worth continuing,
652 ! as there something fundamentally wrong with the restart files
653 if (ierr /= 0) then
654 call restart%close(wfns_file)
655 call restart%close(occ_file)
656 call profiling_out("RESTART_READ")
657 pop_sub(states_elec_load)
658 return
659 end if
660
661 if (states_are_real(st)) then
662 safe_allocate(dpsi(1:mesh%np))
663 else
664 safe_allocate(zpsi(1:mesh%np))
665 end if
666
667 safe_allocate(restart_file(1:st%d%dim, st%st_start:st%st_end, 1:st%nik))
668 safe_allocate(restart_file_present(1:st%d%dim, st%st_start:st%st_end, 1:st%nik))
669 restart_file_present = .false.
670
671 ! Next we read the list of states from the files.
672 ! Errors in reading the information of a specific state from the files are ignored
673 ! at this point, because later we will skip reading the wavefunction of that state.
674 do
675 call restart%read(wfns_file, lines, 1, err)
676 if (err == 0) then
677 read(lines(1), '(a)') char
678 if (char == '%') then
679 !We reached the end of the file
680 exit
681 else
682 read(lines(1), *) ik, char, ist, char, idim, char, filename
683 end if
684 end if
685
686 if (err /= 0 .or. index_is_wrong()) then
687 call restart%read(occ_file, lines, 1, err)
688 cycle
689 end if
690
691 if (ist >= st%st_start .and. ist <= st%st_end .and. &
692 ik >= st%d%kpt%start .and. ik <= st%d%kpt%end) then
693
694 restart_file(idim, ist, ik) = trim(filename)
695 restart_file_present(idim, ist, ik) = .true.
696 end if
697
698 call restart%read(occ_file, lines, 1, err)
699 if (.not. present(lr)) then ! do not read eigenvalues or occupations when reading linear response
700 ! # occupations | eigenvalue[a.u.] | Im(eigenvalue) [a.u.] | k-points | k-weights | filename | ik | ist | idim
701
702 if (err == 0) then
703 read(lines(1), *) my_occ, char, st%eigenval(ist, ik), char, imev, char, &
704 (read_kpoint(idir), char, idir = 1, space%dim), my_kweight
705 ! we do not want to read the k-weights, we have already set them appropriately
706 else
707 ! There is a problem with this states information, so we skip it.
708 if (ist >= st%st_start .and. ist <= st%st_end .and. &
709 ik >= st%d%kpt%start .and. ik <= st%d%kpt%end) then
710 restart_file_present(idim, ist, ik) = .false.
711 end if
712 cycle
713 end if
714
715 kpoint(1:space%dim) = &
716 kpoints%get_point(st%d%get_kpoint_index(ik), absolute_coordinates = .true.)
717 ! FIXME: maybe should ignore ik and just try to match actual vector k-points?
718 if (any(abs(kpoint(1:space%dim) - read_kpoint(1:space%dim)) > 1e-12_real64)) then
719 ! write only once for each k-point so as not to be too verbose
720 if (ist == 1) then
721 write(message(1),'(a,i6)') 'Incompatible restart information: k-point mismatch for ik ', ik
722 write(message(2),'(a,99f18.12)') ' Expected : ', kpoint(1:space%dim)
723 write(message(3),'(a,99f18.12)') ' Read : ', read_kpoint(1:space%dim)
724 call messages_warning(3, namespace=namespace)
725 end if
726 if (ist >= st%st_start .and. ist <= st%st_end .and. &
727 ik >= st%d%kpt%start .and. ik <= st%d%kpt%end) then
728 restart_file_present(idim, ist, ik) = .false.
729 end if
730 end if
731
732 if (read_occ) then
733 st%occ(ist, ik) = my_occ
734 integral_occs = integral_occs .and. &
735 abs((st%occ(ist, ik) - st%smear%el_per_state) * st%occ(ist, ik)) <= m_epsilon
736 end if
737 end if
738 end do
739
740 if (present(iter)) then
741 call restart%read(wfns_file, lines, 1, err)
742 if (err /= 0) then
743 ierr = ierr - 2**8
744 else
745 read(lines(1), *) filename, filename, iter
746 end if
747 end if
748
749 call restart%close(wfns_file)
750 call restart%close(occ_file)
751
752 ! At this point, the restart file should be available
753 ! if not, we fall back to the obf format
754 if (restart_file_format == option__restartfileformatstates__adios2) then
755 if (.not.loct_dir_exists(trim(restart%dir())//"/"//"wavefunctions.bp")) then
756 message(1) = "ADIOS2 restart file not found, falling back to obf format."
757 call messages_warning(1, namespace=namespace)
758 restart_file_format = option__restartfileformatstates__obf
759 end if
760 end if
761
762 ! Now we read the wavefunctions. At this point we need to have all the information from the
763 ! states, occs, and wfns files in order to avoid serialisation of reads, as restart_read
764 ! works as a barrier.
765 safe_allocate(filled(1:st%d%dim, st%st_start:st%st_end, st%d%kpt%start:st%d%kpt%end))
766 filled = .false.
767
768 if (present(lowest_missing)) lowest_missing = st%nst + 1
769
770 iread = 0
771 if (mpi_world%is_root() .and. verbose_) then
772 idone = 0
773 ntodo = st%lnst*st%d%kpt%nlocal*st%d%dim
774 call loct_progress_bar(-1, ntodo)
775 end if
776 if (restart_file_format == option__restartfileformatstates__obf) then
777 do ik = st%d%kpt%start, st%d%kpt%end
778 do ist = st%st_start, st%st_end
779 if (present(skip)) then
780 if (skip(ist)) cycle
781 end if
782
783 do idim = 1, st%d%dim
784
785 if (.not. restart_file_present(idim, ist, ik)) then
786 if (present(lowest_missing)) then
787 lowest_missing(idim, ik) = min(lowest_missing(idim, ik), ist)
788 end if
789 cycle
790 end if
791
792 if (states_are_real(st)) then
793 call restart%read_mesh_function(restart_file(idim, ist, ik), mesh, dpsi, err)
794 if (.not. present(lr)) then
795 call states_elec_set_state(st, mesh, idim, ist, ik, dpsi)
796 else
797 call lalg_copy(mesh%np, dpsi, lr%ddl_psi(:, idim, ist, ik))
798 end if
799 else
800 call restart%read_mesh_function(restart_file(idim, ist, ik), mesh, zpsi, err)
801 if (.not. present(lr)) then
802 call states_elec_set_state(st, mesh, idim, ist, ik, zpsi)
803 else
804 call lalg_copy(mesh%np, zpsi, lr%zdl_psi(:, idim, ist, ik))
805 end if
806 end if
807
808 if (err == 0) then
809 filled(idim, ist, ik) = .true.
810 iread = iread + 1
811 else if (present(lowest_missing)) then
812 lowest_missing(idim, ik) = min(lowest_missing(idim, ik), ist)
813 end if
814
815 if (mpi_world%is_root() .and. verbose_) then
816 idone = idone + 1
817 call loct_progress_bar(idone, ntodo)
818 end if
819
820 end do
821 end do
822 end do
823 else if (restart_file_format == option__restartfileformatstates__adios2) then
824 if (present(lr)) then
825 ! this should never occur because we force obf format for linear response
826 call messages_not_implemented("ADIOS2 restart file format with linear response.")
827 end if
828 if (restart%has_map()) then
829 message(1) = "ADIOS2 restart file format does not support restarting from a different grid."
830 message(2) = "Please set RestartFileFormatStates = obf."
831 call messages_fatal(2, namespace=namespace)
832 end if
833 call states_elec_load_adios2(restart, st, mesh, iread, filled, restart_file_present, ierr, skip, lowest_missing)
834 end if
835
836 safe_deallocate_a(dpsi)
837 safe_deallocate_a(zpsi)
838 safe_deallocate_a(zpsil)
839 safe_deallocate_a(restart_file)
840 safe_deallocate_a(restart_file_present)
841
842 if (mpi_world%is_root() .and. verbose_) then
843 call messages_new_line()
844 end if
845
846 if (st%parallel_in_states .or. st%d%kpt%parallel) then
847 iread_tmp = iread
848 call st%st_kpt_mpi_grp%allreduce(iread_tmp, iread, 1, mpi_integer, mpi_sum)
849 end if
850
851 if (st%d%kpt%parallel) then
852 ! make sure all tasks know lowest_missing over all k-points
853 if (present(lowest_missing)) then
854 safe_allocate(lowest_missing_tmp(1:st%d%dim, 1:st%nik))
855 lowest_missing_tmp = lowest_missing
856 call st%d%kpt%mpi_grp%allreduce(lowest_missing_tmp(1,1), lowest_missing(1,1), st%d%dim*st%nik, &
857 mpi_integer, mpi_min)
858 safe_deallocate_a(lowest_missing_tmp)
859 end if
860 end if
861
862 if (st%restart_fixed_occ .and. iread == st%nst * st%nik * st%d%dim) then
863 ! reset to overwrite whatever smearing may have been set earlier
864 ! only do this if all states could be loaded from restart files
865 call smear_init(st%smear, namespace, st%d%ispin, fixed_occ = .true., integral_occs = integral_occs, kpoints = kpoints)
866 end if
867
868 if (.not. present(lr) .and. .not. present(skip)) call fill_random()
869 ! it is better to initialize lr wfns to zero
870
871 safe_deallocate_a(filled)
872
873 if (ierr == 0 .and. iread /= st%nst * st%nik * st%d%dim) then
874 if (iread > 0) then
875 ierr = iread
876 else
877 ierr = -1
878 end if
879 ! otherwise ierr = 0 would mean either all was read correctly, or nothing at all was read!
880
881 if (.not. present(lr)) then
882 write(str, '(a,i5)') 'Reading states.'
883 else
884 write(str, '(a,i5)') 'Reading states information for linear response.'
885 end if
886 call messages_print_with_emphasis(msg=trim(str), namespace=namespace)
887 if (.not. present(skip)) then
888 write(message(1),'(a,i6,a,i6,a)') 'Only ', iread,' files out of ', &
889 st%nst * st%nik * st%d%dim, ' could be read.'
890 else
891 write(message(1),'(a,i6,a,i6,a)') 'Only ', iread,' files out of ', &
892 st%nst * st%nik * st%d%dim, ' were loaded.'
893 ierr = 0
894 end if
895 call messages_info(1, namespace=namespace)
896 call messages_print_with_emphasis(namespace=namespace)
897 end if
898
899 message(1) = 'Info: States reading done.'
900 if (verbose_) call print_date(trim(message(1))//' ')
901
902 call profiling_out("RESTART_READ")
903 pop_sub(states_elec_load)
904
905 contains
906
907 ! ---------------------------------------------------------
908 subroutine fill_random() !< put random function in orbitals that could not be read.
910
911 do ik = st%d%kpt%start, st%d%kpt%end
912
913 do ist = st%st_start, st%st_end
914 do idim = 1, st%d%dim
915 if (filled(idim, ist, ik)) cycle
916
917 call states_elec_generate_random(st, mesh, kpoints, ist, ist, ik, ik)
918 end do
919 end do
920 end do
921
923 end subroutine fill_random
924
925 ! ---------------------------------------------------------
926
927 logical function index_is_wrong() !< .true. if the index (idim, ist, ik) is not present in st structure...
929
930 if (idim > st%d%dim .or. idim < 1 .or. &
931 ist > st%nst .or. ist < 1 .or. &
932 ik > st%nik .or. ik < 1) then
934 else
935 index_is_wrong = .false.
936 end if
937
939 end function index_is_wrong
940
941 end subroutine states_elec_load
942
943 ! ---------------------------------------------------------
944 subroutine states_elec_load_adios2(restart, st, mesh, iread, filled, restart_file_present, ierr, skip, lowest_missing)
945 type(restart_t), intent(in) :: restart
946 type(states_elec_t), target, intent(inout) :: st
947 class(mesh_t), intent(in) :: mesh
948 integer, intent(out) :: iread
949 logical, intent(out) :: filled(1:, st%st_start:, st%d%kpt%start:)
950 logical, intent(in) :: restart_file_present(:, :, :)
951 integer, intent(out) :: ierr
952 logical, optional, intent(in) :: skip(:)
953 integer, optional, intent(inout) :: lowest_missing(:, :)
954
955#ifdef HAVE_ADIOS2
956 type(adios2_adios) :: adios
957 type(adios2_io) :: io
958 type(adios2_variable) :: var, var_indices
959 type(adios2_attribute) :: attribute
960 type(adios2_engine) :: engine
961 integer :: type_file, type_requested, ib, ik, ndims, pardomains
962 integer :: nst_max
963 integer(int64), allocatable :: var_shape(:)
964 integer(int64), allocatable :: global_indices(:)
965#endif
966
968#ifdef HAVE_ADIOS2
969#ifdef HAVE_MPI
970 call adios2_init(adios, restart%mpi_grp%comm%MPI_VAL, ierr)
971#else
972 call adios2_init(adios, ierr)
973#endif
974 call check_error(.not. adios%valid .or. ierr /= adios2_error_none, "Problem initializing ADIOS2.")
975 call adios2_declare_io(io, adios, "reader", ierr)
976 call check_error(.not. io%valid .or. ierr /= adios2_error_none, "Problem initializing ADIOS2.")
977 ! open restart file
978 call adios2_open(engine, io, trim(restart%dir())//"/"//"wavefunctions.bp", adios2_mode_read, ierr)
979 call check_error(.not. engine%valid .or. ierr /= adios2_error_none, "Problem opening ADIOS2 restart file.")
980
981 call adios2_begin_step(engine, ierr)
982 ! get ParDomains from attribute
983 call adios2_inquire_attribute(attribute, io, "ParDomains", ierr)
984 call check_error(.not. attribute%valid .or. ierr /= adios2_error_none, "Problem inquiring ADIOS2 attribute.")
985 call adios2_attribute_data(pardomains, attribute, ierr)
986 if (states_are_real(st)) then
987 type_requested = adios2_type_dp
988 else
989 type_requested = adios2_type_complex_dp
990 end if
991 call adios2_inquire_variable(var, io, 'wavefunctions', ierr)
992 call check_error(.not. var%valid .or. ierr /= adios2_error_none, "Problem loading ADIOS2 variable (wavefunctions).")
993 ! get variable type and shape
994 call adios2_variable_type(type_file, var, ierr)
995 call adios2_variable_ndims(ndims, var, ierr)
996 call adios2_variable_shape(var_shape, ndims, var, ierr)
997 assert(ndims == 3)
998 nst_max = st%nst
999 if (var_shape(1) < int(st%nst * st%d%dim, int64)) then
1000 ! more states requested than available in the restart data, only read up to nst_max
1001 nst_max = int(var_shape(1), int32) / st%d%dim
1002 end if
1003 if (var_shape(2) /= int(mesh%np_global, int64)) then
1004 ! this should never happen because we check above if there is a map which is only allocated if the grid is different
1005 message(1) = "Error: trying to restart with a different number of grid points. Not supported with ADIOS2 format."
1006 call messages_fatal(1)
1007 end if
1008 if (present(skip)) then
1009 ! read state by state if skip parameter present, otherwise batch by batch which is much better
1010 ! for the performance
1011 call states_elec_adios2_read_state_by_state(st, mesh, engine, var, filled, restart_file_present, iread, skip, lowest_missing)
1012 else
1013 call states_elec_adios2_read_batch_by_batch(st, mesh, engine, var, var_shape, type_requested, type_file, &
1014 nst_max, filled, iread)
1015 end if
1016 ! communicate points for mesh parallelization
1017 if (mesh%mpi_grp%size > 1 .or. pardomains > 1) then
1018 call adios2_inquire_variable(var_indices, io, 'global_indices', ierr)
1019 call check_error(.not. var_indices%valid .or. ierr /= adios2_error_none, &
1020 "Problem loading ADIOS2 variable (global_indices).")
1021 safe_allocate(global_indices(mesh%np))
1022 call adios2_set_selection(var_indices, 1, &
1023 [int(mesh%pv%xlocal-1, int64)], &
1024 [int(mesh%np, int64)], ierr)
1025 call check_error(ierr /= adios2_error_none, "Problem setting ADIOS2 selection (global_indices).")
1026 call adios2_get(engine, var_indices, global_indices, ierr)
1027 call check_error(ierr /= adios2_error_none, "Problem in ADIOS2 get (global_indices).")
1028 end if
1029 call adios2_end_step(engine, ierr)
1030 call check_error(ierr /= adios2_error_none, "Problem in ADIOS2 end_step.")
1031 if (mesh%mpi_grp%size > 1 .or. pardomains > 1) then
1032 do ik = st%d%kpt%start, st%d%kpt%end
1033 do ib = st%group%block_start, st%group%block_end
1034 if (states_are_real(st)) then
1035 call dmesh_batch_exchange_points(mesh, st%group%psib(ib, ik), forward_map=global_indices)
1036 else
1037 call zmesh_batch_exchange_points(mesh, st%group%psib(ib, ik), forward_map=global_indices)
1038 end if
1039 end do
1040 end do
1041 safe_deallocate_a(global_indices)
1042 end if
1043
1044 ! close restart file
1045 call adios2_close(engine, ierr)
1046 call check_error(ierr /= adios2_error_none, "Problem closing ADIOS2 engine.")
1047 call adios2_finalize(adios, ierr)
1048 call check_error(ierr /= adios2_error_none, "Problem finalizing ADIOS2.")
1049#else
1050 ! avoid an empty routine
1051 ierr = 0
1052#endif
1054 end subroutine states_elec_load_adios2
1055
1056
1057#ifdef HAVE_ADIOS2
1058 ! ---------------------------------------------------------
1060 subroutine states_elec_adios2_read_state_by_state(st, mesh, engine, var, filled, restart_file_present, &
1061 iread, skip, lowest_missing)
1062 type(states_elec_t), target, intent(inout) :: st
1063 class(mesh_t), intent(in) :: mesh
1064 type(adios2_engine), intent(inout) :: engine
1065 type(adios2_variable), intent(inout) :: var
1066 logical, intent(inout) :: filled(1:, st%st_start:, st%d%kpt%start:)
1067 logical, intent(in) :: restart_file_present(:, :, :)
1068 integer, intent(inout) :: iread
1069 logical, optional, intent(in) :: skip(:)
1070 integer, optional, intent(inout) :: lowest_missing(:, :)
1071
1072 integer :: ik, ist, idim, ierr
1073 real(real64), allocatable :: dpsi(:)
1074 complex(real64), allocatable :: zpsi(:)
1075
1076 push_sub(states_elec_adios2_read_state_by_state)
1077
1078 do ik = st%d%kpt%start, st%d%kpt%end
1079 if (states_are_real(st)) then
1080 safe_allocate(dpsi(1:mesh%np))
1081 else
1082 safe_allocate(zpsi(1:mesh%np))
1083 end if
1084 do ist = st%st_start, st%st_end
1085 if (present(skip)) then
1086 if (skip(ist)) cycle
1087 end if
1088
1089 do idim = 1, st%d%dim
1090 if (.not. restart_file_present(idim, ist, ik)) then
1091 if (present(lowest_missing)) then
1092 lowest_missing(idim, ik) = min(lowest_missing(idim, ik), ist)
1093 end if
1094 cycle
1095 end if
1096
1097 ! read only one state
1098 call adios2_set_selection(var, 3, &
1099 [int((ist-1)*st%d%dim, int64), int(mesh%pv%xlocal-1, int64), int(ik-1, int64)], &
1100 [1_int64, int(mesh%np, int64), 1_int64], ierr)
1101 ! use sync mode here to immediately read into dpsi/zpsi
1102 if (states_are_real(st)) then
1103 call adios2_get(engine, var, dpsi, adios2_mode_sync, ierr)
1104 call check_error(ierr /= adios2_error_none, "Problem in ADIOS2 get dpsi.")
1105 call states_elec_set_state(st, mesh, idim, ist, ik, dpsi)
1106 else
1107 call adios2_get(engine, var, zpsi, adios2_mode_sync, ierr)
1108 call check_error(ierr /= adios2_error_none, "Problem in ADIOS2 get zpsi.")
1109 call states_elec_set_state(st, mesh, idim, ist, ik, zpsi)
1110 end if
1111
1112 filled(idim, ist, ik) = .true.
1113 iread = iread + 1
1114 end do
1115 end do
1116 safe_deallocate_a(dpsi)
1117 safe_deallocate_a(zpsi)
1118 end do
1119
1120 pop_sub(states_elec_adios2_read_state_by_state)
1121 end subroutine states_elec_adios2_read_state_by_state
1123 ! ---------------------------------------------------------
1125 subroutine states_elec_adios2_read_batch_by_batch(st, mesh, engine, var, var_shape, type_requested, type_file, &
1126 nst_max, filled, iread)
1127 type(states_elec_t), target, intent(inout) :: st
1128 class(mesh_t), intent(in) :: mesh
1129 type(adios2_engine), intent(inout) :: engine
1130 type(adios2_variable), intent(inout) :: var
1131 integer(int64), intent(in) :: var_shape(:)
1132 integer, intent(in) :: type_requested
1133 integer, intent(in) :: type_file
1134 integer, intent(in) :: nst_max
1135 logical, intent(inout) :: filled(1:, st%st_start:, st%d%kpt%start:)
1136 integer, intent(inout) :: iread
1137
1138 integer :: ib, ik, ist, ip, adios2_mode, nst_linear, ierr
1139 class(wfs_elec_t), pointer :: psib
1140 real(real64), allocatable :: dff_pack(:, :)
1141
1142 push_sub(states_elec_adios2_read_batch_by_batch)
1143
1144 do ik = st%d%kpt%start, st%d%kpt%end
1145 ! if there are more k points requested than available, skip the ones that are not available
1146 ! this can, e.g., happen for unocc calculations with a k points path
1147 if (ik > var_shape(3)) cycle
1148 do ib = st%group%block_start, st%group%block_end
1149 select case (st%group%psib(ib, ik)%status())
1150 case (batch_not_packed)
1151 safe_allocate(psib)
1152 call st%group%psib(ib, ik)%copy_to(psib, copy_data=.false.)
1153 call psib%do_pack(batch_packed, copy=.false.)
1154 adios2_mode = adios2_mode_sync
1155 case (batch_packed)
1156 psib => st%group%psib(ib, ik)
1157 !adios2_mode = adios2_mode_deferred
1158 ! switch to sync mode for the time being, seems to solve some errors
1159 ! at the moment, it is unclear why
1160 adios2_mode = adios2_mode_sync
1161 case (batch_device_packed)
1162 safe_allocate(psib)
1163 call st%group%psib(ib, ik)%copy_to(psib)
1164 call psib%do_unpack(force=.true., copy=.false.)
1165 ! make sure we have a packed batch
1166 call psib%do_pack(batch_packed)
1167 adios2_mode = adios2_mode_sync
1168 end select
1169 assert(psib%status() == batch_packed)
1170
1171 ! read only up to the highest available state in the restart data
1172 nst_linear = st%group%psib(ib, ik)%nst_linear
1173 if (states_elec_block_max(st, ib) > nst_max) then
1174 nst_linear = (nst_max - states_elec_block_min(st, ib) + 1) * st%d%dim
1175 end if
1176 if (nst_linear <= 0) cycle
1177
1178 ! select part of variable to write, correspond to the memory of the current batch
1179 ! first argument gives the start in the global array
1180 ! second argument gives the size of the local batch
1181 call adios2_set_selection(var, 3, &
1182 [int((states_elec_block_min(st, ib)-1)*st%d%dim, int64), int(mesh%pv%xlocal-1, int64), int(ik-1, int64)], &
1183 [int(nst_linear, int64), int(mesh%np, int64), 1_int64], ierr)
1184 call check_error(ierr /= adios2_error_none, "Problem setting ADIOS2 selection.")
1185
1186 if (type_requested == adios2_type_dp .and. type_file == adios2_type_dp) then
1187 ! set memory layout of one batch, take into account padding and ghost/boundary points
1188 call adios2_set_memory_selection(var, 3, [0_int64, 0_int64, 0_int64], &
1189 [int(size(psib%dff_pack, 1), int64), &
1190 int(size(psib%dff_pack, 2), int64), 1_int64], ierr)
1191 call check_error(ierr /= adios2_error_none, "Problem setting ADIOS2 memory selection (reading real).")
1192
1193 call adios2_get(engine, var, psib%dff_pack, adios2_mode, ierr)
1194 call check_error(ierr /= adios2_error_none, "Problem in ADIOS2 get (reading real).")
1195 else if (type_requested == adios2_type_complex_dp .and. type_file == adios2_type_complex_dp) then
1196 ! set memory layout of one batch, take into account padding and ghost/boundary points
1197 call adios2_set_memory_selection(var, 3, [0_int64, 0_int64, 0_int64], &
1198 [int(size(psib%zff_pack, 1), int64), &
1199 int(size(psib%zff_pack, 2), int64), 1_int64], ierr)
1200 call check_error(ierr /= adios2_error_none, "Problem setting ADIOS2 memory selection (reading complex).")
1201
1202 call adios2_get(engine, var, psib%zff_pack, adios2_mode, ierr)
1203 call check_error(ierr /= adios2_error_none, "Problem in ADIOS2 get (reading complex).")
1204 else if (type_requested == adios2_type_complex_dp .and. type_file == adios2_type_dp) then
1205 ! read into a local real variable and copy to the complex batch buffer
1206 safe_allocate(dff_pack(size(psib%zff_pack, 1), size(psib%zff_pack, 2)))
1207 ! set memory layout of one batch, take into account padding and ghost/boundary points
1208 call adios2_set_memory_selection(var, 3, [0_int64, 0_int64, 0_int64], &
1209 [int(size(psib%zff_pack, 1), int64), &
1210 int(size(psib%zff_pack, 2), int64), 1_int64], ierr)
1211 call check_error(ierr /= adios2_error_none, "Problem setting ADIOS2 memory selection (reading complex from real).")
1212
1213 ! use sync mode here to immediately read into dff_pack
1214 call adios2_get(engine, var, dff_pack, adios2_mode_sync, ierr)
1215 call check_error(ierr /= adios2_error_none, "Problem in ADIOS2 get (reading complex from real).")
1216 ! copy to batch memory
1217 !$omp parallel do private(ist)
1218 do ip = 1, mesh%np
1219 do ist = 1, nst_linear
1220 psib%zff_pack(ist, ip) = cmplx(dff_pack(ist, ip), m_zero, real64)
1221 end do
1222 end do
1223 safe_deallocate_a(dff_pack)
1224 else
1225 message(1) = "Error: requested to read complex states for restarting, but calculation has real states."
1226 call messages_fatal(1)
1227 end if
1228 ! make sure iread and filled are set correctly
1229 iread = iread + nst_linear
1230 do ist = 1, nst_linear
1231 filled(st%group%psib(ib, ik)%linear_to_idim(ist), st%group%psib(ib, ik)%linear_to_ist(ist), ik) = .true.
1232 end do
1233
1234 select case (st%group%psib(ib, ik)%status())
1235 case (batch_not_packed)
1236 call psib%do_unpack()
1237 call psib%copy_data_to(mesh%np, st%group%psib(ib, ik))
1238 call psib%end()
1239 safe_deallocate_p(psib)
1240 case (batch_packed)
1241 nullify(psib)
1242 case (batch_device_packed)
1243 call psib%do_pack()
1244 call psib%copy_data_to(mesh%np, st%group%psib(ib, ik))
1245 call psib%end()
1246 safe_deallocate_p(psib)
1247 end select
1248 end do
1249 end do
1250
1251 pop_sub(states_elec_adios2_read_batch_by_batch)
1252 end subroutine states_elec_adios2_read_batch_by_batch
1253
1254 ! ---------------------------------------------------------
1255 ! only used for adios2 error checking
1256 subroutine check_error(condition, error_message)
1257 logical, intent(in) :: condition
1258 character(len=*), intent(in) :: error_message
1259 if (condition) then
1260 message(1) = error_message
1261 call messages_fatal(1)
1262 end if
1263 end subroutine check_error
1264#endif
1265
1266 subroutine states_elec_dump_rho(restart, st, mesh, ierr, iter)
1267 type(restart_t), intent(in) :: restart
1268 type(states_elec_t), intent(in) :: st
1269 class(mesh_t), intent(in) :: mesh
1270 integer, intent(out) :: ierr
1271 integer, optional, intent(in) :: iter
1272
1273 integer :: iunit, isp, err, err2(2)
1274 character(len=MAX_PATH_LEN) :: filename
1275 character(len=300) :: lines(2)
1276
1277 push_sub(states_elec_dump_rho)
1278
1279 ierr = 0
1280
1281 if (restart%skip()) then
1282 pop_sub(states_elec_dump_rho)
1283 return
1284 end if
1285
1286 message(1) = "Debug: Writing density restart."
1287 call messages_info(1, debug_only=.true.)
1288
1290
1291 !write the densities
1292 iunit = restart%open('density')
1293 lines(1) = '# #spin #nspin filename'
1294 lines(2) = '%densities'
1295 call restart%write(iunit, lines, 2, err)
1296 if (err /= 0) ierr = ierr + 1
1297
1298 err2 = 0
1299 do isp = 1, st%d%nspin
1300 filename = get_filename_with_spin('density', st%d%nspin, isp)
1301 write(lines(1), '(i8,a,i8,a)') isp, ' | ', st%d%nspin, ' | "'//trim(adjustl(filename))//'"'
1302 call restart%write(iunit, lines, 1, err)
1303 if (err /= 0) err2(1) = err2(1) + 1
1304
1305 call restart%write_mesh_function(filename, mesh, st%rho(:,isp), err)
1306 if (err /= 0) err2(2) = err2(2) + 1
1307
1308 end do
1309 if (err2(1) /= 0) ierr = ierr + 2
1310 if (err2(2) /= 0) ierr = ierr + 4
1311
1312 lines(1) = '%'
1313 call restart%write(iunit, lines, 1, err)
1314 if (err /= 0) ierr = ierr + 8
1315 if (present(iter)) then
1316 write(lines(1),'(a,i7)') 'Iter = ', iter
1317 call restart%write(iunit, lines, 1, err)
1318 if (err /= 0) ierr = ierr + 16
1319 end if
1320 call restart%close(iunit)
1321
1323
1324 message(1) = "Debug: Writing density restart done."
1325 call messages_info(1, debug_only=.true.)
1326
1327 pop_sub(states_elec_dump_rho)
1328 end subroutine states_elec_dump_rho
1329
1330
1331 ! ---------------------------------------------------------
1332 subroutine states_elec_load_rho(restart, st, mesh, ierr)
1333 type(restart_t), intent(in) :: restart
1334 type(states_elec_t), intent(inout) :: st
1335 class(mesh_t), intent(in) :: mesh
1336 integer, intent(out) :: ierr
1337
1338 integer :: err, err2, isp
1339 character(len=MAX_PATH_LEN) :: filename
1340
1341 push_sub(states_elec_load_rho)
1342
1343 ierr = 0
1344
1345 if (restart%skip()) then
1346 ierr = -1
1347 pop_sub(states_elec_load_rho)
1348 return
1349 end if
1350
1351 message(1) = "Debug: Reading density restart."
1352 call messages_info(1, debug_only=.true.)
1353
1354 ! skip for now, since we know what the files are going to be called
1355 !read the densities
1356! iunit_rho = io_open(trim(dir)//'/density', restart%namespace, action='write')
1357! call iopar_read(st%dom_st_kpt_mpi_grp, iunit_rho, line, err)
1358! call iopar_read(st%dom_st_kpt_mpi_grp, iunit_rho, line, err)
1359! we could read the iteration 'iter' too, not sure if that is useful.
1360
1361 err2 = 0
1362 do isp = 1, st%d%nspin
1363 filename = get_filename_with_spin('density', st%d%nspin, isp)
1364! if (st%dom_st_kpt_mpi_grp%is_root()) then
1365! read(iunit_rho, '(i8,a,i8,a)') isp, ' | ', st%d%nspin, ' | "'//trim(adjustl(filename))//'"'
1366! end if
1367 call restart%read_mesh_function(filename, mesh, st%rho(:,isp), err)
1368 if (err /= 0) err2 = err2 + 1
1369
1370 end do
1371 if (err2 /= 0) ierr = ierr + 1
1372
1373 message(1) = "Debug: Reading density restart done."
1374 call messages_info(1, debug_only=.true.)
1375
1376 pop_sub(states_elec_load_rho)
1377 end subroutine states_elec_load_rho
1378
1379 subroutine states_elec_dump_frozen(restart, space, st, mesh, ierr)
1380 type(restart_t), intent(in) :: restart
1381 class(space_t), intent(in) :: space
1382 type(states_elec_t), intent(in) :: st
1383 class(mesh_t), intent(in) :: mesh
1384 integer, intent(out) :: ierr
1385
1386 integer :: isp, err, err2(2), idir
1387 character(len=MAX_PATH_LEN) :: filename
1388
1389 push_sub(states_elec_dump_frozen)
1390
1391 ierr = 0
1392
1393 assert(allocated(st%frozen_rho))
1394
1395 if (restart%skip()) then
1397 return
1398 end if
1399
1400 message(1) = "Debug: Writing frozen densities restart."
1401 call messages_info(1, debug_only=.true.)
1402
1404
1405 err2 = 0
1406 do isp = 1, st%d%nspin
1407 filename = get_filename_with_spin('frozen_rho', st%d%nspin, isp)
1408
1409 call restart%write_mesh_function(filename, mesh, st%frozen_rho(:,isp), err)
1410 if (err /= 0) err2(2) = err2(2) + 1
1411
1412 if (allocated(st%frozen_tau)) then
1413 filename = get_filename_with_spin('frozen_tau', st%d%nspin, isp)
1414 call restart%write_mesh_function(filename, mesh, st%frozen_tau(:,isp), err)
1415 if (err /= 0) err2 = err2 + 1
1416 end if
1417
1418 if (allocated(st%frozen_gdens)) then
1419 do idir = 1, space%dim
1420 if (st%d%nspin == 1) then
1421 write(filename, fmt='(a,i1)') 'frozen_gdens-dir', idir
1422 else
1423 write(filename, fmt='(a,i1,a,i1)') 'frozen_tau-dir', idir, '-', isp
1424 end if
1425 call restart%write_mesh_function(filename, mesh, st%frozen_gdens(:,idir,isp), err)
1426 if (err /= 0) err2 = err2 + 1
1427 end do
1428 end if
1429
1430 if (allocated(st%frozen_ldens)) then
1431 filename = get_filename_with_spin('frozen_ldens', st%d%nspin, isp)
1432 call restart%write_mesh_function(filename, mesh, st%frozen_ldens(:,isp), err)
1433 if (err /= 0) err2 = err2 + 1
1434 end if
1435
1436 end do
1437 if (err2(1) /= 0) ierr = ierr + 2
1438 if (err2(2) /= 0) ierr = ierr + 4
1439
1441
1442 message(1) = "Debug: Writing frozen densities restart done."
1443 call messages_info(1, debug_only=.true.)
1444
1446 end subroutine states_elec_dump_frozen
1447
1448
1449 ! ---------------------------------------------------------
1450 subroutine states_elec_load_frozen(restart, space, st, mesh, ierr)
1451 type(restart_t), intent(in) :: restart
1452 class(space_t), intent(in) :: space
1453 type(states_elec_t), intent(inout) :: st
1454 class(mesh_t), intent(in) :: mesh
1455 integer, intent(out) :: ierr
1456
1457 integer :: err, err2, isp, idir
1458 character(len=MAX_PATH_LEN) :: filename
1459
1460 push_sub(states_elec_load_frozen)
1461
1462 assert(allocated(st%frozen_rho))
1463
1464 ierr = 0
1465
1466 if (restart%skip()) then
1467 ierr = -1
1469 return
1470 end if
1471
1472 message(1) = "Debug: Reading densities restart."
1473 call messages_info(1, debug_only=.true.)
1474
1475 err2 = 0
1476 do isp = 1, st%d%nspin
1477 filename = get_filename_with_spin('frozen_rho', st%d%nspin, isp)
1478 call restart%read_mesh_function(filename, mesh, st%frozen_rho(:,isp), err)
1479 if (err /= 0) err2 = err2 + 1
1480
1481 if (allocated(st%frozen_tau)) then
1482 filename = get_filename_with_spin('frozen_tau', st%d%nspin, isp)
1483 call restart%read_mesh_function(filename, mesh, st%frozen_tau(:,isp), err)
1484 if (err /= 0) err2 = err2 + 1
1485 end if
1486
1487 if (allocated(st%frozen_gdens)) then
1488 do idir = 1, space%dim
1489 if (st%d%nspin == 1) then
1490 write(filename, fmt='(a,i1)') 'frozen_gdens-dir', idir
1491 else
1492 write(filename, fmt='(a,i1,a,i1)') 'frozen_tau-dir', idir, '-', isp
1493 end if
1494 call restart%read_mesh_function(filename, mesh, st%frozen_gdens(:,idir,isp), err)
1495 if (err /= 0) err2 = err2 + 1
1496 end do
1497 end if
1498
1499 if (allocated(st%frozen_ldens)) then
1500 filename = get_filename_with_spin('frozen_ldens', st%d%nspin, isp)
1501 call restart%read_mesh_function(filename, mesh, st%frozen_ldens(:,isp), err)
1502 if (err /= 0) err2 = err2 + 1
1503 end if
1504
1505 end do
1506 if (err2 /= 0) ierr = ierr + 1
1507
1508 message(1) = "Debug: Reading frozen densities restart done."
1509 call messages_info(1, debug_only=.true.)
1510
1512 end subroutine states_elec_load_frozen
1513
1514
1515 ! ---------------------------------------------------------
1518 subroutine states_elec_read_user_def_orbitals(mesh, namespace, space, st)
1519 class(mesh_t), intent(in) :: mesh
1520 type(namespace_t), intent(in) :: namespace
1521 class(space_t), intent(in) :: space
1522 type(states_elec_t), intent(inout) :: st
1523
1524 type(block_t) :: blk
1525 integer :: ip, id, is, ik, nstates, state_from, ierr, ncols
1526 integer :: ib, idim, inst, inik, normalize
1527 real(real64) :: xx(space%dim), rr, psi_re, psi_im
1528 character(len=150) :: filename
1529 complex(real64), allocatable :: zpsi(:, :)
1530
1531 integer, parameter :: &
1532 state_from_formula = 1, &
1533 state_from_file = -10010, &
1534 normalize_yes = 1, &
1535 normalize_no = 0
1536
1538
1539 !%Variable UserDefinedStates
1540 !%Type block
1541 !%Section States
1542 !%Description
1543 !% Instead of using the ground state as initial state for
1544 !% time-propagations it might be interesting in some cases
1545 !% to specify alternate states. Like with user-defined
1546 !% potentials, this block allows you to specify formulas for
1547 !% the orbitals at <i>t</i>=0.
1548 !%
1549 !% Example:
1550 !%
1551 !% <tt>%UserDefinedStates
1552 !% <br>&nbsp;&nbsp; 1 | 1 | 1 | formula | "exp(-r^2)*exp(-i*0.2*x)" | normalize_yes
1553 !% <br>%</tt>
1554 !%
1555 !% The first column specifies the component of the spinor,
1556 !% the second column the number of the state and the third
1557 !% contains <i>k</i>-point and spin quantum numbers. Column four
1558 !% indicates that column five should be interpreted as a formula
1559 !% for the corresponding orbital.
1560 !%
1561 !% Alternatively, if column four states <tt>file</tt> the state will
1562 !% be read from the file given in column five.
1563 !%
1564 !% <tt>%UserDefinedStates
1565 !% <br>&nbsp;&nbsp; 1 | 1 | 1 | file | "/path/to/file" | normalize_no
1566 !% <br>%</tt>
1567 !%
1568 !% Octopus reads first the ground-state orbitals from
1569 !% the <tt>restart/gs</tt> directory. Only the states that are
1570 !% specified in the above block will be overwritten with
1571 !% the given analytic expression for the orbital.
1572 !%
1573 !% The sixth (optional) column indicates whether <tt>Octopus</tt> should renormalize
1574 !% the orbital. The default (no sixth column given) is to renormalize.
1575 !%
1576 !%Option file -10010
1577 !% Read initial orbital from file.
1578 !% Accepted file formats, detected by extension: obf, ncdf and csv (real only).
1579 !% For csv files, the following formatting rules apply:
1580 !%
1581 !% The functions in this file read an array from an ascii matrix (csv) file.
1582 !% Format with values "valueXYZ" as follows
1583 !%
1584 !% File values.csv:
1585 !% --------
1586 !% value111 value211 value311 value411
1587 !% value121 value221 value321 value421
1588 !% value131 value231 value331 value431
1589 !%
1590 !% value112 value212 value312 value412
1591 !% value122 value222 value322 value422
1592 !% value132 value232 value332 value432
1593 !%
1594 !% value113 value213 value313 value413
1595 !% value123 value223 value323 value423
1596 !% value133 value233 value333 value433
1597 !% --------
1598 !%
1599 !% That is, every line scans the x-coordinate, every XY-plane as a table of values
1600 !% and all XY-planes separated by an empty row.
1601 !%
1602 !% The given matrix is interpolated/stretched to fit the calculation
1603 !% box defined in input file.
1604 !%
1605 !% Note that there must not be any empty line at the end of the file.
1606 !%
1607 !% Calculation box shape must be "parallelepiped".
1608 !%
1609 !% The delimiter can be a tab, a comma or a space.
1610 !%Option formula 1
1611 !% Calculate initial orbital by given analytic expression.
1612 !%Option normalize_yes 1
1613 !% Normalize orbitals (default).
1614 !%Option normalize_no 0
1615 !% Do not normalize orbitals.
1616 !%End
1617 if (parse_block(namespace, 'UserDefinedStates', blk) == 0) then
1618
1619 call messages_print_with_emphasis(msg=trim('Substitution of orbitals'), namespace=namespace)
1620
1621 ! find out how many lines (i.e. states) the block has
1622 nstates = parse_block_n(blk)
1623
1624 safe_allocate(zpsi(1:mesh%np, 1:st%d%dim))
1625
1626 ! read all lines
1627 do ib = 1, nstates
1628 ! Check that number of columns is five or six.
1629 ncols = parse_block_cols(blk, ib - 1)
1630 if (ncols < 5 .or. ncols > 6) then
1631 message(1) = 'Each line in the UserDefinedStates block must have'
1632 message(2) = 'five or six columns.'
1633 call messages_fatal(2, namespace=namespace)
1634 end if
1635
1636 call parse_block_integer(blk, ib - 1, 0, idim)
1637 call parse_block_integer(blk, ib - 1, 1, inst)
1638 call parse_block_integer(blk, ib - 1, 2, inik)
1639
1640 ! Calculate from expression or read from file?
1641 call parse_block_integer(blk, ib - 1, 3, state_from)
1642
1643 ! loop over all states
1644 do id = 1, st%d%dim
1645 do is = 1, st%nst
1646 do ik = 1, st%nik
1647
1648 if (.not.(id == idim .and. is == inst .and. ik == inik )) cycle
1649
1650 ! We first parse the value, such that the parsing and stdout looks correct
1651 ! However, if rank 0 is not responsible for the parsing, the variables in the formula
1652 ! will still appear as not parsed (pure k-point/state parallelization case).
1653 ! This could be solved by a dummy call to parse_expression
1654 select case (state_from)
1655 case (state_from_formula)
1656 ! parse formula string
1657 call parse_block_string( &
1658 blk, ib - 1, 4, st%user_def_states(id, is, ik))
1659
1660 write(message(1), '(a,3i5)') 'Substituting state of orbital with k, ist, dim = ', ik, is, id
1661 write(message(2), '(2a)') ' with the expression:'
1662 write(message(3), '(2a)') ' ',trim(st%user_def_states(id, is, ik))
1663 call messages_info(3, namespace=namespace)
1664
1665 ! convert to C string
1666 call conv_to_c_string(st%user_def_states(id, is, ik))
1667
1668 case (state_from_file)
1669 ! The input format can be coded in column four now. As it is
1670 ! not used now, we just say "file".
1671 ! Read the filename.
1672 call parse_block_string(blk, ib - 1, 4, filename)
1673
1674 write(message(1), '(a,3i5)') 'Substituting state of orbital with k, ist, dim = ', ik, is, id
1675 write(message(2), '(2a)') ' with data from file:'
1676 write(message(3), '(2a)') ' ',trim(filename)
1677 call messages_info(3, namespace=namespace)
1678 end select
1679
1680
1681 ! does the block entry match and is this node responsible?
1682 if (.not.(st%st_start <= is .and. st%st_end >= is &
1683 .and. st%d%kpt%start <= ik .and. st%d%kpt%end >= ik)) cycle
1684
1685 select case (state_from)
1686
1687 case (state_from_formula)
1688 ! fill states with user-defined formulas
1689 do ip = 1, mesh%np
1690 xx = mesh%x(ip, :)
1691 rr = norm2(xx)
1692
1693 ! parse user-defined expressions
1694 call parse_expression(psi_re, psi_im, space%dim, xx, rr, m_zero, st%user_def_states(id, is, ik))
1695 ! fill state
1696 zpsi(ip, 1) = cmplx(psi_re, psi_im, real64)
1697 end do
1698
1699 case (state_from_file)
1700 ! finally read the state
1701 call zio_function_input(filename, namespace, space, mesh, zpsi(:, 1), ierr)
1702 if (ierr > 0) then
1703 message(1) = 'Could not read the file!'
1704 write(message(2),'(a,i1)') 'Error code: ', ierr
1705 call messages_fatal(2, namespace=namespace)
1706 end if
1707
1708 case default
1709 message(1) = 'Wrong entry in UserDefinedStates, column 4.'
1710 message(2) = 'You may state "formula" or "file" here.'
1711 call messages_fatal(2, namespace=namespace)
1712 end select
1713
1714 ! normalize orbital
1715 if (parse_block_cols(blk, ib - 1) == 6) then
1716 call parse_block_integer(blk, ib - 1, 5, normalize)
1717 else
1718 normalize = 1
1719 end if
1720 select case (normalize)
1721 case (normalize_no)
1722 call states_elec_set_state(st, mesh, id, is, ik, zpsi(:, 1))
1723 case (normalize_yes)
1724 assert(st%d%dim == 1)
1725 call zmf_normalize(mesh, st%d%dim, zpsi)
1726 call states_elec_set_state(st, mesh, is, ik, zpsi)
1727 case default
1728 message(1) = 'The sixth column in UserDefinedStates may either be'
1729 message(2) = '"normalize_yes" or "normalize_no"'
1730 call messages_fatal(2, namespace=namespace)
1731 end select
1732
1733 end do
1734 end do
1735 end do
1736
1737 end do
1738
1739 safe_deallocate_a(zpsi)
1740
1741 call parse_block_end(blk)
1742 call messages_print_with_emphasis(namespace=namespace)
1743
1744 else
1745 call messages_variable_is_block(namespace, 'UserDefinedStates')
1746 end if
1747
1750
1751 ! ---------------------------------------------------------
1752 ! This is needed as for the generalized Bloch theorem we need to label
1753 ! the states from the expectation value of Sz computed from the GS.
1754 subroutine states_elec_dump_spin(restart, st, ierr)
1755 type(restart_t), intent(in) :: restart
1756 type(states_elec_t), intent(in) :: st
1757 integer, intent(out) :: ierr
1758
1759 integer :: iunit_spin
1760 integer :: err, err2(2), ik, ist
1761 character(len=300) :: lines(3)
1762
1763 push_sub(states_elec_dump_spin)
1764
1765 ierr = 0
1766
1767 if (restart%skip()) then
1768 pop_sub(states_elec_dump_spin)
1769 return
1770 end if
1771
1772 call profiling_in("RESTART_WRITE_SPIN")
1773
1775
1776 iunit_spin = restart%open('spin')
1777 lines(1) = '# #k-point #st #spin(x) spin(y) spin(z)'
1778 call restart%write(iunit_spin, lines, 1, err)
1779 if (err /= 0) ierr = ierr + 1
1780
1781 err2 = 0
1782 do ik = 1, st%nik
1783 do ist = 1, st%nst
1784 write(lines(1), '(i8,a,i8,3(a,f18.12))') ik, ' | ', ist, ' | ', &
1785 st%spin(1,ist,ik), ' | ', st%spin(2,ist,ik),' | ', st%spin(3,ist,ik)
1786 call restart%write(iunit_spin, lines, 1, err)
1787 if (err /= 0) err2(1) = err2(1) + 1
1788 end do ! st%nst
1789 end do ! st%nik
1790 lines(1) = '%'
1791 call restart%write(iunit_spin, lines, 1, err)
1792 if (err2(1) /= 0) ierr = ierr + 8
1793 if (err2(2) /= 0) ierr = ierr + 16
1794
1795 call restart%close(iunit_spin)
1796
1798
1799 call profiling_out("RESTART_WRITE_SPIN")
1800 pop_sub(states_elec_dump_spin)
1801 end subroutine states_elec_dump_spin
1802
1803
1804
1805
1806 ! ---------------------------------------------------------
1811 subroutine states_elec_load_spin(restart, st, ierr)
1812 type(restart_t), intent(in) :: restart
1813 type(states_elec_t), intent(inout) :: st
1814 integer, intent(out) :: ierr
1815
1816 integer :: spin_file, err, ik, ist
1817 character(len=256) :: lines(3)
1818 real(real64) :: spin(3)
1819 character(len=1) :: char
1820
1821
1822 push_sub(states_elec_load_spin)
1823
1824 ierr = 0
1825
1826 if (restart%skip()) then
1827 ierr = -1
1828 pop_sub(states_elec_load_spin)
1829 return
1830 end if
1831
1832 call profiling_in("RESTART_READ_SPIN")
1833
1834 ! open files to read
1835 spin_file = restart%open('spin')
1836 ! Skip two lines.
1837 call restart%read(spin_file, lines, 1, err)
1838 if (err /= 0) ierr = ierr - 2**7
1839
1840 ! If any error occured up to this point then it is not worth continuing,
1841 ! as there something fundamentally wrong with the restart files
1842 if (ierr /= 0) then
1843 call restart%close(spin_file)
1844 call profiling_out("RESTART_READ_SPIN")
1845 pop_sub(states_elec_load_spin)
1846 return
1847 end if
1848
1849 ! Next we read the list of states from the files.
1850 ! Errors in reading the information of a specific state from the files are ignored
1851 ! at this point, because later we will skip reading the wavefunction of that state.
1852 do
1853 call restart%read(spin_file, lines, 1, err)
1854 read(lines(1), '(a)') char
1855 if (char == '%') then
1856 !We reached the end of the file
1857 exit
1858 else
1859 read(lines(1), *) ik, char, ist, char, spin(1), char, spin(2), char, spin(3)
1860 end if
1861
1862 if (err /= 0) cycle
1863
1864 st%spin(1:3, ist, ik) = spin(1:3)
1865 end do
1866
1867 call restart%close(spin_file)
1868
1869 call profiling_out("RESTART_READ_SPIN")
1870 pop_sub(states_elec_load_spin)
1871 end subroutine states_elec_load_spin
1872
1873 ! ---------------------------------------------------------
1874 subroutine states_elec_transform(st, namespace, space, restart, mesh, kpoints, prefix)
1875 type(states_elec_t), intent(inout) :: st
1876 type(namespace_t), intent(in) :: namespace
1877 class(space_t), intent(in) :: space
1878 type(restart_t), intent(inout) :: restart
1879 class(mesh_t), intent(in) :: mesh
1880 type(kpoints_t), intent(in) :: kpoints
1881 character(len=*), optional, intent(in) :: prefix
1882
1883 type(states_elec_t) :: stin
1884 type(block_t) :: blk
1885 complex(real64), allocatable :: rotation_matrix(:,:), psi(:, :)
1886 integer :: ist, jst, ncols, iqn
1887 character(len=256) :: block_name
1888
1889 push_sub(states_elec_transform)
1890
1891 block_name = trim(optional_default(prefix, "")) // "TransformStates"
1892
1893 !%Variable TransformStates
1894 !%Type block
1895 !%Default no
1896 !%Section States
1897 !%Description
1898 !% Before starting the <tt>td</tt> calculation, the initial states (that are
1899 !% read from the <tt>restart/gs</tt> directory, which should have been
1900 !% generated in a previous ground-state calculation) can be "transformed"
1901 !% among themselves. The block <tt>TransformStates</tt> gives the transformation matrix
1902 !% to be used. The number of rows and columns of the matrix should equal the number
1903 !% of the states present in the time-dependent calculation (the independent
1904 !% spin and <i>k</i>-point subspaces are all transformed equally); the number of
1905 !% columns should be equal to the number of states present in the
1906 !% <tt>restart/gs</tt> directory. This number may be different: for example,
1907 !% one could have run previously in <tt>unocc</tt> mode in order to obtain unoccupied
1908 !% Kohn-Sham states, and therefore <tt>restart/gs</tt> will contain more states.
1909 !% These states can be used in the transformation.
1910 !%
1911 !% Note that the code will not check the orthonormality of the new states!
1912 !%
1913 !% Each line provides the coefficients of the new states, in terms of
1914 !% the old ones. The coefficients are complex, but the imaginary part will be
1915 !% ignored for real wavefunctions.
1916 !% Note: This variable cannot be used when parallel in states.
1917 !%End
1918 if (parse_is_defined(namespace, trim(block_name))) then
1919 if (parse_block(namespace, trim(block_name), blk) == 0) then
1920 if (st%parallel_in_states) then
1921 call messages_not_implemented(trim(block_name) // " parallel in states", namespace=namespace)
1922 end if
1923 if (parse_block_n(blk) /= st%nst) then
1924 message(1) = "Number of rows in block " // trim(block_name) // " must equal number of states in this calculation."
1925 call messages_fatal(1, namespace=namespace)
1926 end if
1927 call states_elec_copy(stin, st, exclude_wfns = .true.)
1928 call states_elec_look_and_load(restart, namespace, space, stin, mesh, kpoints)
1929
1930 ! FIXME: rotation matrix should be R_TYPE
1931 safe_allocate(rotation_matrix(1:stin%nst, 1:stin%nst))
1932 safe_allocate(psi(1:mesh%np, 1:st%d%dim))
1933
1934 rotation_matrix = diagonal_matrix(stin%nst, m_z1)
1935
1936 do ist = 1, st%nst
1937 ncols = parse_block_cols(blk, ist-1)
1938 if (ncols /= stin%nst) then
1939 write(message(1),'(a,i6,a,i6,3a,i6,a)') "Number of columns (", ncols, ") in row ", ist, " of block ", &
1940 trim(block_name), " must equal number of states (", stin%nst, ") read from gs restart."
1941 call messages_fatal(1, namespace=namespace)
1942 end if
1943 do jst = 1, stin%nst
1944 call parse_block_cmplx(blk, ist - 1, jst - 1, rotation_matrix(jst, ist))
1945 end do
1946 end do
1947
1948 call parse_block_end(blk)
1949
1950 do iqn = st%d%kpt%start, st%d%kpt%end
1951 if (states_are_real(st)) then
1952 call states_elec_rotate(stin, mesh, real(rotation_matrix, real64) , iqn)
1953 else
1954 call states_elec_rotate(stin, mesh, rotation_matrix, iqn)
1955 end if
1956
1957 do ist = st%st_start, st%st_end
1958 call states_elec_get_state(stin, mesh, ist, iqn, psi)
1959 call states_elec_set_state(st, mesh, ist, iqn, psi)
1960 end do
1961
1962 end do
1963
1964 safe_deallocate_a(rotation_matrix)
1965 safe_deallocate_a(psi)
1966
1967 call states_elec_end(stin)
1968
1969 else
1970 call messages_variable_is_block(namespace, trim(block_name))
1971 end if
1972
1973 end if
1974
1975 pop_sub(states_elec_transform)
1976 end subroutine states_elec_transform
1977
1979
1980
1981!! Local Variables:
1982!! mode: f90
1983!! coding: utf-8
1984!! End:
ssize_t read(int __fd, void *__buf, size_t __nbytes) __attribute__((__access__(__write_only__
long int random(void)
Definition: getopt_f.c:1407
Copies a vector x, to a vector y.
Definition: lalg_basic.F90:188
block signals while writing the restart files
Definition: restart.F90:318
unblock signals when writing restart is finished
Definition: restart.F90:325
This module implements batches of mesh functions.
Definition: batch.F90:135
integer, parameter, public batch_not_packed
functions are stored in CPU memory, unpacked order
Definition: batch.F90:286
integer, parameter, public batch_device_packed
functions are stored in device memory in packed order
Definition: batch.F90:286
integer, parameter, public batch_packed
functions are stored in CPU memory, in transposed (packed) order
Definition: batch.F90:286
This module handles the calculation mode.
integer, parameter, public p_strategy_domains
parallelization in domains
integer, parameter, public spinors
real(real64), parameter, public m_zero
Definition: global.F90:191
real(real64), parameter, public m_epsilon
Definition: global.F90:207
complex(real64), parameter, public m_z1
Definition: global.F90:202
subroutine, public zio_function_input(filename, namespace, space, mesh, ff, ierr, map)
Reads a mesh function from file filename, and puts it into ff. If the map argument is passed,...
System information (time, memory, sysname)
Definition: loct.F90:117
subroutine, public loct_progress_bar(a, maxcount)
A wrapper around the progress bar, such that it can be silenced without needing to dress the call wit...
Definition: loct.F90:276
logical function, public loct_dir_exists(dirname)
Definition: loct.F90:349
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:117
This module defines functions over batches of mesh functions.
Definition: mesh_batch.F90:118
subroutine, public zmesh_batch_exchange_points(mesh, aa, forward_map, backward_map)
This functions exchanges points of a mesh according to a certain map. Two possible maps can be given....
subroutine, public dmesh_batch_exchange_points(mesh, aa, forward_map, backward_map)
This functions exchanges points of a mesh according to a certain map. Two possible maps can be given....
This module defines various routines, operating on mesh functions.
subroutine, public zmf_normalize(mesh, dim, psi, norm)
Normalize a mesh function psi.
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:120
subroutine, public messages_print_with_emphasis(msg, iunit, namespace)
Definition: messages.F90:898
subroutine, public messages_not_implemented(feature, namespace)
Definition: messages.F90:1091
character(len=512), private msg
Definition: messages.F90:167
subroutine, public messages_variable_is_block(namespace, name)
Definition: messages.F90:1047
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:525
subroutine, public print_date(str)
Definition: messages.F90:983
subroutine, public messages_new_line()
Definition: messages.F90:1112
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:162
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:410
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:594
type(mpi_grp_t), public mpi_world
Definition: mpi.F90:272
This module handles the communicators for the various parallelization strategies.
Definition: multicomm.F90:147
this module contains the low-level part of the output system
Definition: output_low.F90:117
character(len=max_path_len) function, public get_filename_with_spin(output, nspin, spin_index)
Returns the filame as output, or output-spX is spin polarized.
Definition: output_low.F90:233
Some general things and nomenclature:
Definition: par_vec.F90:173
integer(int64) function, public par_vec_local2global(pv, ip)
Returns global index of local point ip.
Definition: par_vec.F90:805
logical function, public parse_is_defined(namespace, name)
Definition: parser.F90:455
subroutine, public parse_block_string(blk, l, c, res, convert_to_c)
Definition: parser.F90:810
integer function, public parse_block(namespace, name, blk, check_varinfo_)
Definition: parser.F90:615
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
Definition: profiling.F90:631
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
Definition: profiling.F90:554
subroutine, public smear_init(this, namespace, ispin, fixed_occ, integral_occs, kpoints)
Definition: smear.F90:187
pure logical function, public states_are_complex(st)
pure logical function, public states_are_real(st)
This module handles spin dimensions of the states and the k-point distribution.
integer pure function, public states_elec_block_max(st, ib)
return index of last state in block ib
subroutine, public states_elec_end(st)
finalize the states_elec_t object
subroutine, public states_elec_allocate_wfns(st, mesh, wfs_type, skip, packed)
Allocates the KS wavefunctions defined within a states_elec_t structure.
subroutine, public states_elec_copy(stout, stin, exclude_wfns, exclude_eigenval, special)
make a (selective) copy of a states_elec_t object
subroutine, public states_elec_generate_random(st, mesh, kpoints, ist_start_, ist_end_, ikpt_start_, ikpt_end_, normalized)
randomize states
integer pure function, public states_elec_block_min(st, ib)
return index of first state in block ib
subroutine, public states_elec_look(restart, nik, dim, nst, ierr)
Reads the 'states' file in the restart directory, and finds out the nik, dim, and nst contained in it...
This module handles reading and writing restart information for the states_elec_t.
subroutine, public states_elec_read_user_def_orbitals(mesh, namespace, space, st)
the routine reads formulas for user-defined wavefunctions from the input file and fills the respectiv...
subroutine, public states_elec_load_frozen(restart, space, st, mesh, ierr)
subroutine, public states_elec_look_and_load(restart, namespace, space, st, mesh, kpoints, is_complex, packed)
subroutine states_elec_dump_adios2(restart, st, mesh, ierr)
subroutine, public states_elec_transform(st, namespace, space, restart, mesh, kpoints, prefix)
subroutine, public states_elec_dump(restart, space, st, mesh, kpoints, ierr, iter, lr, verbose)
subroutine, public states_elec_load(restart, namespace, space, st, mesh, kpoints, ierr, iter, lr, lowest_missing, label, verbose, skip)
returns in ierr: <0 => Fatal error, or nothing read =0 => read all wavefunctions >0 => could only rea...
subroutine, public states_elec_load_spin(restart, st, ierr)
returns in ierr: <0 => Fatal error, or nothing read =0 => read all wavefunctions >0 => could only rea...
subroutine, public states_elec_dump_frozen(restart, space, st, mesh, ierr)
subroutine states_elec_load_adios2(restart, st, mesh, iread, filled, restart_file_present, ierr, skip, lowest_missing)
subroutine, public states_elec_load_rho(restart, st, mesh, ierr)
subroutine, public states_elec_dump_rho(restart, st, mesh, ierr, iter)
subroutine, public states_elec_dump_spin(restart, st, ierr)
subroutine, public conv_to_c_string(str)
converts to c string
Definition: string.F90:257
type(type_t), public type_float
Definition: types.F90:135
type(type_t), public type_cmplx
Definition: types.F90:136
subroutine fill_random()
logical function index_is_wrong()
Describes mesh distribution to nodes.
Definition: mesh.F90:187
The states_elec_t class contains all electronic wave functions.
batches of electronic states
Definition: wfs_elec.F90:141
int true(void)