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