Octopus
convert.F90
Go to the documentation of this file.
1!! Copyright (C) 2013 J. Alberdi-Rodriguez, J. Jornet-Somoza
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
20
21#include "global.h"
22
23program oct_convert
24 use batch_oct_m
27 use debug_oct_m
29 use fft_oct_m
31 use global_oct_m
32 use io_oct_m
35 use ions_oct_m
36 use kick_oct_m
38 use loct_oct_m
40 use mesh_oct_m
41 use mpi_oct_m
44 use output_oct_m
46 use parser_oct_m
50 use space_oct_m
52 use string_oct_m
53 use unit_oct_m
55 use utils_oct_m
56
57 implicit none
58
59 character(len=256) :: config_str
60 integer :: ierr
61
62 call getopt_init(ierr)
63 config_str = trim(get_config_opts()) // trim(get_optional_libraries())
64 if (ierr == 0) call getopt_octopus(config_str)
65 call getopt_end()
66
67 call global_init()
68
69 call parser_init()
70
71 call messages_init()
72
73 call io_init()
75 call messages_experimental("oct-convert utility")
76
77 call print_header()
78
79 call messages_print_with_emphasis(msg="Convert mode", namespace=global_namespace)
81
84
85 call convert()
86
87 call fft_all_end()
89 call io_end()
90 call print_date("Calculation ended on ")
91 call messages_end()
92
93 call parser_end()
94
95 call global_end()
96
97contains
98
99 ! -------------
104 subroutine convert()
105 type(electrons_t), pointer :: sys
106
107 character(MAX_PATH_LEN) :: basename, folder, ref_name, ref_folder, folder_default
108 integer :: c_start, c_end, c_step, c_start_default, length, c_how
109 logical :: iterate_folder, subtract_file
110 integer, parameter :: CONVERT_FORMAT = 1, fourier_transform = 2, operation = 3
111
112 push_sub(convert)
113
114 call calc_mode_par%set_parallelization(p_strategy_states, default = .false.)
115 sys => electrons_t(global_namespace, generate_epot=.false.)
116 call sys%init_parallelization(mpi_world)
117
118 message(1) = 'Info: Converting files'
119 message(2) = ''
120 call messages_info(2)
121
122 !%Variable ConvertFilename
123 !%Type string
124 !%Default "density"
125 !%Section Utilities::oct-convert
126 !%Description
127 !% Input filename. The original filename which is going to be converted in the format
128 !% specified in <tt>OutputFormat</tt>. It is going to convert various files, it should
129 !% only contain the beginning of the name. For instance, in the case of the restart
130 !% files, it should be one space ' '.
131 !%End
132 call parse_variable(global_namespace, 'ConvertFilename', 'density', basename)
133 if (basename == " ") basename = ""
134 ! Delete the extension if present
135 length = len_trim(basename)
136 if (length > 4) then
137 if (basename(length-3:length) == '.obf') then
138 basename = trim(basename(1:length-4))
139 end if
140 end if
141
142 !%Variable ConvertHow
143 !%Type integer
144 !%Default convert_format
145 !%Section Utilities::oct-convert
146 !%Description
147 !% Select how the mesh function will be converted.
148 !%Option format 1
149 !% The format of the mesh function will be convert from the binary file.obf.
150 !% The format of the output function is set by OutputHow variable.
151 !%Option fourier_transform 2
152 !% A fourier transform of the mesh function will be computed.
153 !% It requieres that ConvertStart and ConvertEnd have to be set.
154 !%Option operation 3
155 !% Convert utility will generate a new mesh function constructed by linear
156 !% combination of scalar function of different mesh functions,
157 !%End
158 call parse_variable(global_namespace, 'ConvertHow', convert_format, c_how)
159
160 !%Variable ConvertIterateFolder
161 !%Type logical
162 !%Default true
163 !%Section Utilities::oct-convert
164 !%Description
165 !% This variable decides if a folder is going to be iterated or the
166 !% filename is going to be iterated.
167 !%End
168 call parse_variable(global_namespace, 'ConvertIterateFolder', .true., iterate_folder)
169
170 if (iterate_folder) then
171 folder_default = 'output_iter/td.'
172 c_start_default = 0
173 else
174 folder_default = 'restart'
175 c_start_default = 1
176 end if
177
178 !%Variable ConvertFolder
179 !%Type string
180 !%Section Utilities::oct-convert
181 !%Description
182 !% The folder name where the input files are. The default is
183 !% <tt>output_iter/td.</tt> if <tt>ConvertIterateFolder = true</tt>, otherwise <tt>restart</tt>.
184 !%End
185 call parse_variable(global_namespace, 'ConvertFolder', folder_default, folder)
186 call add_last_slash(folder)
187
188 !%Variable ConvertStart
189 !%Type integer
190 !%Section Utilities::oct-convert
191 !%Description
192 !% The starting number of the filename or folder.
193 !% Default is 0 if <tt>ConvertIterateFolder = true</tt>, otherwise 1.
194 !%End
195 call parse_variable(global_namespace, 'ConvertStart', c_start_default, c_start)
196
197 !%Variable ConvertEnd
198 !%Type integer
199 !%Default 1
200 !%Section Utilities::oct-convert
201 !%Description
202 !% The last number of the filename or folder.
203 !%End
204 call parse_variable(global_namespace, 'ConvertEnd', 1, c_end)
205
206 !%Variable ConvertStep
207 !%Type integer
208 !%Default 1
209 !%Section Utilities::oct-convert
210 !%Description
211 !% The padding between the filenames or folder.
212 !%End
213 call parse_variable(global_namespace, 'ConvertStep', 1, c_step)
214
215 !%Variable ConvertSubtractFilename
216 !%Type string
217 !%Default density
218 !%Section Utilities::oct-convert
219 !%Description
220 !% Input filename. The file which is going to subtracted to rest of the files.
221 !%End
222 call parse_variable(global_namespace, 'ConvertSubtractFilename', 'density', ref_name)
223 if (ref_name == " ") ref_name = ""
224 ! Delete the extension if present
225 length = len_trim(ref_name)
226 if (length > 4) then
227 if (ref_name(length-3:length) == '.obf') then
228 ref_name = trim(ref_name(1:length-4))
229 end if
230 end if
231
232 !%Variable ConvertSubtract
233 !%Type logical
234 !%Default false
235 !%Section Utilities::oct-convert
236 !%Description
237 !% Decides if a reference file is going to be subtracted.
238 !%End
239 call parse_variable(global_namespace, 'ConvertSubtract', .false., subtract_file)
240
241 !%Variable ConvertSubtractFolder
242 !%Type string
243 !%Default .
244 !%Section Utilities::oct-convert
245 !%Description
246 !% The folder name which is going to be subtracted.
247 !%End
248 call parse_variable(global_namespace, 'ConvertSubtractFolder', '.', ref_folder)
249 call add_last_slash(folder)
250
251 select case (c_how)
252 CASE(operation)
253 call convert_operate(sys%gr, global_namespace, sys%space, sys%ions, sys%mc, sys%outp)
254
255 CASE(fourier_transform)
256 ! Compute Fourier transform
257 call convert_transform(sys%gr, global_namespace, sys%space, sys%ions, sys%mc, sys%kpoints, basename, folder, &
258 c_start, c_end, c_step, sys%outp, subtract_file, &
259 ref_name, ref_folder)
260
261 CASE(convert_format)
262 call convert_low(sys%gr, global_namespace, sys%space, sys%ions, sys%hm%psolver, sys%mc, basename, folder, &
263 c_start, c_end, c_step, sys%outp, iterate_folder, &
264 subtract_file, ref_name, ref_folder)
265 end select
266
267 safe_deallocate_p(sys)
268
269 pop_sub(convert)
270 end subroutine convert
271
272 ! ---------------------------------------------------------
275 subroutine convert_low(mesh, namespace, space, ions, psolver, mc, basename, in_folder, c_start, c_end, c_step, outp, &
276 iterate_folder, subtract_file, ref_name, ref_folder)
277 class(mesh_t), intent(in) :: mesh
278 type(namespace_t), intent(in) :: namespace
279 class(space_t), intent(in) :: space
280 type(ions_t), intent(in) :: ions
281 type(poisson_t), intent(in) :: psolver
282 type(multicomm_t), intent(in) :: mc
283 character(len=*), intent(inout) :: basename
284 character(len=*), intent(in) :: in_folder
285 integer, intent(in) :: c_start
286 integer, intent(in) :: c_end
287 integer, intent(in) :: c_step
288 type(output_t), intent(in) :: outp
289 logical, intent(in) :: iterate_folder
291 logical, intent(in) :: subtract_file
292 character(len=*), intent(inout) :: ref_name
293 character(len=*), intent(inout) :: ref_folder
294
295 type(restart_t) :: restart
296 integer :: ierr, ii, folder_index, output_i
297 character(MAX_PATH_LEN) :: filename, out_name, folder, frmt, restart_folder
298 real(real64), allocatable :: read_ff(:), read_rff(:), pot(:)
299
300 push_sub(convert_low)
301
302 safe_allocate(read_ff(1:mesh%np))
303 safe_allocate(read_rff(1:mesh%np))
304 safe_allocate(pot(1:mesh%np))
305 read_rff(:) = m_zero
306
307 write(message(1),'(5a,i5,a,i5,a,i5)') "Converting '", trim(in_folder), "/", trim(basename), &
308 "' from ", c_start, " to ", c_end, " every ", c_step
309 call messages_info(1)
310
311 if (subtract_file) then
312 write(message(1),'(a,a,a,a)') "Reading ref-file from ", trim(ref_folder), trim(ref_name),".obf"
313 call restart%init(namespace, restart_undefined, restart_type_load, mc, ierr, &
314 dir=trim(ref_folder), mesh = mesh)
315 ! FIXME: why only real functions? Please generalize.
316 if (ierr == 0) then
317 call restart%read_mesh_function(space, trim(ref_name), mesh, read_rff, ierr)
318 call restart%end()
319 else
320 write(message(1),'(2a)') "Failed to read from ref-file ", trim(ref_name)
321 write(message(2), '(2a)') "from folder ", trim(ref_folder)
322 call messages_fatal(2)
323 end if
324 end if
325
326 ! Initialize the restart directory from <tt>ConvertFolder</tt> value.
327 ! This directory has to have the files 'grid' and 'lxyz.obf'
328 ! and the files that are going to be converged, must be inside this folder
329 if (iterate_folder) then
330 ! Delete the last / and find the previous /, if any
331 folder = in_folder(1:len_trim(in_folder)-1)
332 folder_index = index(folder, '/', .true.)
333 restart_folder = folder(1:folder_index)
334 else
335 restart_folder = in_folder
336 end if
337 call restart%init(namespace, restart_undefined, restart_type_load, mc, ierr, &
338 dir=trim(restart_folder), mesh = mesh)
339 call loct_progress_bar(-1, c_end-c_start)
340 do ii = c_start, c_end, c_step
341 if (iterate_folder) then
342 ! Delete the last / and add the corresponding folder number
343 write(folder,'(a,i0.7,a)') in_folder(folder_index+1:len_trim(in_folder)-1),ii,"/"
344 write(filename, '(a,a,a)') trim(folder), trim(basename)
345 out_name = trim(basename)
346 else
347 folder = ""
348 if (c_start /= c_end) then
349 ! Here, we are only considering 10 character long filenames.
350 ! Subtract the initial part given at 'ConvertFilename' from the format and pad
351 ! with zeros.
352 write(frmt,'(a,i0,a)')"(a,i0.",10-len_trim(basename),")"
353 write(filename, fmt=trim(frmt)) trim(basename), ii
354 write(out_name, '(a)') trim(filename)
355 else
356 ! Assuming filename is given complete in the 'ConvertFilename'
357 write(filename, '(a,a,a,a)') trim(folder),"/", trim(basename)
358 filename = basename
359 write(out_name, '(a)') trim(basename)
360 end if
361 end if
362
363 ! Read the obf file
364 call restart%read_mesh_function(space, trim(filename), mesh, read_ff, ierr)
365
366 if (ierr /= 0) then
367 write(message(1), '(a,a)') "Error reading the file ", trim(filename)
368 write(message(2), '(a,i4)') "Error code: ",ierr
369 write(message(3), '(a)') "Skipping...."
371 cycle
372 end if
373 if (subtract_file) then
374 read_ff(:) = read_ff(:) - read_rff(:)
375 write(out_name, '(a,a)') trim(out_name),"-ref"
376 end if
377 ! Write the corresponding output
378 do output_i = lbound(outp%how, 1), ubound(outp%how, 1)
379 if (outp%how(output_i) /= 0) then
380 call dio_function_output(outp%how(output_i), trim(restart_folder)//trim(folder), &
381 trim(out_name), namespace, space, mesh, read_ff, units_out%length**(-space%dim), ierr, &
382 pos=ions%pos, atoms=ions%atom)
383 end if
384 end do
385 if (outp%what(option__output__potential)) then
386 write(out_name, '(a)') "potential"
387 call dpoisson_solve(psolver, namespace, pot, read_ff)
388 call dio_function_output(outp%how(option__output__potential), trim(restart_folder)//trim(folder), &
389 trim(out_name), namespace, space, mesh, pot, units_out%energy, ierr, pos=ions%pos, atoms=ions%atom)
390 end if
391 call loct_progress_bar(ii-c_start, c_end-c_start)
392 ! It does not matter if the current write has failed for the next iteration
393 ierr = 0
394 end do
395 call restart%end()
396
397 safe_deallocate_a(read_ff)
398 safe_deallocate_a(read_rff)
399 safe_deallocate_a(pot)
400 pop_sub(convert_low)
401 end subroutine convert_low
402
403 ! ---------------------------------------------------------
406 subroutine convert_transform(mesh, namespace, space, ions, mc, kpoints, basename, in_folder, c_start, c_end, c_step, outp, &
407 subtract_file, ref_name, ref_folder)
408 class(mesh_t), intent(in) :: mesh
409 type(namespace_t), intent(in) :: namespace
410 class(space_t), intent(in) :: space
411 type(ions_t), intent(in) :: ions
412 type(multicomm_t), intent(in) :: mc
413 type(kpoints_t), intent(in) :: kpoints
414 character(len=*), intent(inout) :: basename
415 character(len=*), intent(in) :: in_folder
416 integer, intent(in) :: c_start
417 integer, intent(in) :: c_end
418 integer, intent(in) :: c_step
419 type(output_t), intent(in) :: outp
420 logical, intent(in) :: subtract_file
421 character(len=*), intent(inout) :: ref_name
422 character(len=*), intent(inout) :: ref_folder
423
424 integer :: ierr, i_space, i_time, nn(1:3), optimize_parity(1:3), wd_info, output_i
425 integer :: i_energy, e_end, e_start, e_point, chunk_size, read_count, t_point
426 logical :: optimize(1:3)
427 integer :: folder_index
428 character(MAX_PATH_LEN) :: filename, folder, restart_folder
429 real(real64) :: fdefault, w_max
430 real(real64), allocatable :: read_ft(:), read_rff(:), point_tmp(:,:)
431
432 integer, parameter :: FAST_FOURIER = 1, standard_fourier = 2
433 type(kick_t) :: kick
434 integer :: ft_method
435 type(spectrum_t) :: spectrum
436 type(batch_t) :: tdrho_b, wdrho_b
437 real(real64), allocatable :: tdrho_a(:,:,:), wdrho_a(:,:,:)
438 type(fft_t) :: fft
439
440 type(restart_t) :: restart
441 complex(real64), allocatable :: out_fft(:)
442
443 real(real64) :: start_time
444 integer :: time_steps
445 real(real64) :: dt
446 real(real64) :: dw
447 real(real64) :: max_energy
448 real(real64) :: min_energy
449
450 push_sub(convert_transform)
451
452 ! set default time_step as dt from TD section
453 fdefault = m_zero
454 call parse_variable(namespace, 'TDTimeStep', fdefault, dt, unit = units_inp%time)
455 if (dt <= m_zero) then
456 write(message(1),'(a)') 'Input: TDTimeStep must be positive.'
457 write(message(2),'(a)') 'Input: TDTimeStep reset to 0. Check input file'
458 call messages_info(2)
459 end if
460
461 call io_mkdir('wd.general', namespace)
462 wd_info = io_open('wd.general/wd.info', global_namespace, action='write')
463 call messages_print_with_emphasis(msg="Fourier Transform Options", iunit=wd_info)
464
465 !%Variable ConvertEnergyMin
466 !%Type float
467 !%Default 0.0
468 !%Section Utilities::oct-convert
469 !%Description
470 !% Minimum energy to output from Fourier transform.
471 !%End
472 call parse_variable(namespace, 'ConvertEnergyMin', m_zero, min_energy, units_inp%energy)
473 call messages_print_var_value('ConvertEnergyMin', min_energy, unit = units_out%energy, iunit=wd_info)
474
475 !%Variable ConvertReadSize
476 !%Type integer
477 !%Default mesh%np
478 !%Section Utilities::oct-convert
479 !%Description
480 !% How many points are read at once. For the parallel run this has not been
481 !% yet tested, so it should be one. For the serial run, a number
482 !% of 100-1000 will speed-up the execution time by this factor.
483 !%End
484 call parse_variable(namespace, 'ConvertReadSize', mesh%np, chunk_size)
485 call messages_print_var_value('ConvertReadSize', chunk_size, iunit=wd_info)
486 !Check that default value is set when ConvertReadSize = 0
487 if (chunk_size == 0) chunk_size = mesh%np
488 ! Parallel version just work in domains and chunk_size equal to mesh%np
489 if (mesh%mpi_grp%size > 1 .and. chunk_size /= mesh%np) then
490 write(message(1),*)'Incompatible value for ConvertReadSize and Parallelizaion in Domains'
491 write(message(2),*)'Use the default value for ConvertReadSize (or set it to 0)'
492 call messages_fatal(2)
493 end if
494
495 ! Calculate the limits in frequency space.
496 start_time = c_start * dt
497 dt = dt * c_step
498 time_steps = (c_end - c_start) / c_step
499 dw = m_two * m_pi / (dt * time_steps)
500 w_max = m_two * m_pi / dt
502 !%Variable ConvertEnergyMax
503 !%Type float
504 !%Default w_max
505 !%Section Utilities::oct-convert
506 !%Description
507 !% Maximum energy to output from Fourier transform.
508 !%End
509 fdefault = units_from_atomic(units_inp%energy, w_max)
510 call parse_variable(namespace, 'ConvertEnergyMax',fdefault, max_energy, units_inp%energy)
511 if (max_energy > w_max) then
512 write(message(1),'(a,f12.7)')'Impossible to set ConvertEnergyMax to ', &
513 units_from_atomic(units_inp%energy, max_energy)
514 write(message(2),'(a)')'ConvertEnergyMax is too large.'
515 write(message(3),'(a,f12.7,a)')'ConvertEnergyMax reset to ', &
516 units_from_atomic(units_inp%energy, w_max),'[' // trim(units_abbrev(units_out%energy)) // ']'
517 call messages_info(3)
518 max_energy = w_max
519 end if
520 call messages_print_var_value('ConvertEnergyMax', max_energy, unit = units_out%energy, iunit=wd_info)
521
522 !%Variable ConvertFTMethod
523 !%Type integer
524 !%Default FAST_FOURIER
525 !%Section Utilities::oct-convert
526 !%Description
527 !% Describes the method used to perform the Fourier Transform
528 !%Option fast_fourier 1
529 !% Uses Fast Fourier Transform as implemented in the external library.
530 !%Option standard_fourier 2
531 !% Uses polinomial approach to the computation of discrete Fourier Transform.
532 !% It uses the same variable described in how to obtain spectrum from
533 !% a time-propagation calculation.
534 !%End
535 call parse_variable(namespace, 'ConvertFTMethod', 1, ft_method)
536 call messages_print_var_option('ConvertFTMethod', ft_method, iunit=wd_info)
537
538 !TODO: check if e_point can be used instead of e_point+1
539 safe_allocate(read_ft(0:time_steps))
540 read_ft = m_zero
541 safe_allocate(read_rff(1:mesh%np))
542
543 select case (ft_method)
544 case (fast_fourier)
545 nn(1) = time_steps + 1
546 nn(2) = 1
547 nn(3) = 1
548 safe_allocate(out_fft(0:time_steps))
549 optimize = .false.
550 optimize_parity = -1
551 call fft_init(fft, nn, 1, fft_real, fftlib_fftw, optimize, optimize_parity)
552 case (standard_fourier)
553 !%Variable ConvertEnergyStep
554 !%Type float
555 !%Default <math>2 \pi / T</math>, where <math>T</math> is the total propagation time
556 !%Section Utilities::oct-convert
557 !%Description
558 !% Energy step to output from Fourier transform.
559 !% Sampling rate for the Fourier transform. If you supply a number equal or smaller than zero, then
560 !% the sampling rate will be <math>2 \pi / T</math>, where <math>T</math> is the total propagation time.
561 !%End
562 fdefault = m_two * m_pi / (dt * time_steps)
563 call parse_variable(namespace, 'ConvertEnergyStep',fdefault, dw, units_inp%energy)
564 if (dw <= m_zero) dw = m_two * m_pi / (dt * time_steps)
565
566 call spectrum_init(spectrum, namespace, dw, w_max)
567 ! Manually setting already defined variables on spectrum.
568 spectrum%start_time = c_start * dt
569 spectrum%end_time = c_end * dt
570 spectrum%energy_step = dw
571 spectrum%max_energy = max_energy
572 safe_allocate(tdrho_a(0:time_steps, 1, 1))
573 safe_allocate(wdrho_a(0:time_steps, 1, 1))
574 end select
575 call messages_print_var_value('ConvertEnergyStep', dw, unit = units_out%energy, iunit=wd_info)
576
577 !TODO: set system variable common for all the program in
578 ! order to use call kick_init(kick, sy%st%d%nspin, sys%space%dim, sys%ions%periodic_dim)
579 call kick_init(kick, namespace, space, kpoints, 1)
580
581 e_start = nint(min_energy / dw)
582 e_end = nint(max_energy / dw)
583 write(message(1),'(a,1x,i0.7,a,f12.7,a,i0.7,a,f12.7,a)')'Frequency index:',e_start,'(',&
584 units_from_atomic(units_out%energy, e_start * dw),')-',e_end,'(',units_from_atomic(units_out%energy, e_end * dw),')'
585 write(message(2),'(a,f12.7,a)')'Frequency Step, dw: ', units_from_atomic(units_out%energy, dw), &
586 '[' // trim(units_abbrev(units_out%energy)) // ']'
587 call messages_info(2)
588
589 if (subtract_file) then
590 write(message(1),'(a,a,a,a)') "Reading ref-file from ", trim(ref_folder), trim(ref_name),".obf"
591 call messages_info(1)
592 call restart%init(namespace, restart_undefined, restart_type_load, mc, ierr, &
593 dir=trim(ref_folder), mesh = mesh)
594 ! FIXME: why only real functions? Please generalize.
595 if (ierr == 0) then
596 call restart%read_mesh_function(space, trim(ref_name), mesh, read_rff, ierr)
597 call restart%end()
598 else
599 write(message(1),'(2a)') "Failed to read from ref-file ", trim(ref_name)
600 write(message(2), '(2a)') "from folder ", trim(ref_folder)
601 call messages_fatal(2)
602 end if
603 end if
604
605 call messages_print_with_emphasis(msg="File Information", iunit=wd_info)
606 do i_energy = e_start, e_end
607 write(filename,'(a14,i0.7,a1)')'wd.general/wd.',i_energy,'/'
608 write(message(1),'(a,a,f12.7,a,1x,i7,a)')trim(filename),' w =', &
609 units_from_atomic(units_out%energy,(i_energy) * dw), &
610 '[' // trim(units_abbrev(units_out%energy)) // ']'
611 if (mpi_world%is_root()) write(wd_info,'(a)') message(1)
612 call messages_info(1)
613 call io_mkdir(trim(filename), namespace)
614 end do
615 call io_close(wd_info)
616
617 if (mesh%parallel_in_domains) then
618 ! Delete the last / and find the previous /, if any
619 folder = in_folder(1:len_trim(in_folder)-1)
620 folder_index = index(folder, '/', .true.)
621 restart_folder = folder(1:folder_index)
622 call restart%init(namespace, restart_undefined, restart_type_load, mc, ierr, &
623 dir=trim(restart_folder), mesh = mesh)
624 end if
625
626 !For each mesh point, open density file and read corresponding point.
627 if (mpi_world%is_root()) call loct_progress_bar(-1, mesh%np)
628
629 safe_allocate(point_tmp(1:chunk_size, 0:time_steps))
630 read_count = 0
631 ! Space
632 do i_space = 1, mesh%np
633 ! Time
634 t_point = 0
635 do i_time = c_start, c_end, c_step
636
637 if (mesh%parallel_in_domains .and. i_space == 1) then
638 write(folder,'(a,i0.7,a)') in_folder(folder_index+1:len_trim(in_folder)-1),i_time,"/"
639 write(filename, '(a,a,a)') trim(folder), trim(basename)
640 call restart%read_mesh_function(space, trim(filename), mesh, point_tmp(:, t_point), ierr)
641 else
642 ! Here, we always iterate folders
643 ! Delete the last / and add the corresponding folder number
644 write(folder,'(a,i0.7,a)') in_folder(1:len_trim(in_folder)-1),i_time,"/"
645 write(filename, '(a,a,a,a)') trim(folder), trim(basename), ".obf"
646 if (mod(i_space-1, chunk_size) == 0) then
647 call profiling_in("READING")
648 !TODO: check for any error on the whole file before reading by parts.
649 call io_binary_read(trim(filename), i4_to_i8(chunk_size), point_tmp(1:chunk_size, t_point), &
650 ierr, offset=i4_to_i8(i_space-1))
651 call profiling_out("READING")
652 if (i_time == c_start) read_count = 0
653 end if
654 end if
655
656 if (ierr /= 0 .and. i_space == 1) then
657 write(message(1), '(a,a,2i10)') "Error reading the file ", trim(filename), i_space, i_time
658 write(message(2), '(a)') "Skipping...."
659 write(message(3), '(a,i0)') "Error :", ierr
660 call messages_warning(3)
661 cycle
662 end if
663
664 if (i_time == c_start) read_count = read_count + 1
665 if (subtract_file) then
666 read_ft(t_point) = point_tmp(read_count, t_point) - read_rff(i_space)
667 else
668 read_ft(t_point) = point_tmp(read_count, t_point)
669 end if
670
671 t_point = t_point + 1
672 end do ! Time
673
674 select case (ft_method)
675 case (fast_fourier)
676 call profiling_in("CONVERT_FFTW")
677 call dfft_forward(fft, read_ft, out_fft)
678 call profiling_out("CONVERT_FFTW")
679 ! Should the value be multiplied by dt ??? as in standard discrete Fourier Transform ?
680 point_tmp(read_count, 0:time_steps) = aimag(out_fft(0:time_steps)) * dt
681 case (standard_fourier)
682 tdrho_a(0:time_steps, 1, 1) = read_ft(0:time_steps)
683 call batch_init(tdrho_b, 1, 1, 1, tdrho_a)
684 call batch_init(wdrho_b, 1, 1, 1, wdrho_a)
685 call spectrum_signal_damp(spectrum%damp, spectrum%damp_factor, c_start + 1, c_start + time_steps + 1, &
686 kick%time, dt, tdrho_b)
687 call spectrum_fourier_transform(spectrum%method, spectrum%transform, spectrum%noise, &
688 c_start + 1, c_start + time_steps + 1, kick%time, dt, tdrho_b, min_energy, max_energy, &
689 spectrum%energy_step, wdrho_b)
690 call tdrho_b%end()
691 call wdrho_b%end()
692 do e_point = e_start, e_end
693 point_tmp(read_count, e_point) = - wdrho_a(e_point, 1, 1)
694 end do
695 end select
696
697 if (mod(i_space-1, 1000) == 0 .and. mpi_world%is_root()) then
698 call loct_progress_bar(i_space-1, mesh%np)
699 end if
700
701 !print out wd densities from (ii-chunksize,ii] if running in serial
702 if (mesh%mpi_grp%size == 1) then
703 if (mod(i_space, chunk_size) == 0) then
704 write(message(1),'(a)') ""
705 write(message(2),'(a,i0)') "Writing binary output: step ", i_space/chunk_size
706 call messages_info(2)
707 do i_energy = e_start, e_end
708 write(filename,'(a14,i0.7,a12)')'wd.general/wd.',i_energy,'/density.obf'
709 ! If it is the first time entering here, write the header. But, only once
710 if (i_space == chunk_size) then
711 call dwrite_header(trim(filename), mesh%np_global, ierr)
712 end if
713 call io_binary_write(trim(filename), i4_to_i8(chunk_size), point_tmp(1:chunk_size, i_energy), ierr, &
714 nohead = .true.)
715 end do
716 end if
717 end if
718 end do ! Space
719
720 if (mpi_world%is_root()) call loct_progress_bar(mesh%np, mesh%np)
721 call mesh%mpi_grp%barrier()
722
723 if (mesh%parallel_in_domains) then
724 do i_energy = e_start, e_end
725 write(filename,'(a14,i0.7,a1)')'wd.general/wd.',i_energy,'/'
726 do output_i = lbound(outp%how, 1), ubound(outp%how, 1)
727 if (outp%how(output_i) /= 0) then
728 call dio_function_output(0_8, trim(filename), &
729 trim('density'), namespace, space, mesh, point_tmp(:, i_energy), &
730 units_out%length**(-space%dim), ierr, pos=ions%pos, atoms=ions%atom)
731 end if
732 end do
733 end do
734 call restart%end()
735 else
736 ! write the output files
737 if (any(outp%how /= option__outputformat__binary)) then
738 do i_energy = e_start, e_end
739 write(filename,'(a14,i0.7,a1)')'wd.general/wd.',i_energy,'/'
740 call io_binary_read(trim(filename)//'density.obf', i4_to_i8(mesh%np), read_rff, ierr)
741 do output_i = lbound(outp%how, 1), ubound(outp%how, 1)
742 if ((outp%how(output_i) /= 0) .and. (outp%how(output_i) /= option__outputformat__binary)) then
743 call dio_function_output(outp%how(output_i), trim(filename), &
744 trim('density'), namespace, space, mesh, read_rff, &
745 units_out%length**(-space%dim), ierr, pos=ions%pos, atoms=ions%atom)
746 end if
747 end do
748 end do
749 end if
750 end if
751
752 safe_deallocate_a(point_tmp)
753 safe_deallocate_a(read_ft)
754 safe_deallocate_a(read_rff)
755
756 select case (ft_method)
757 case (fast_fourier)
758 safe_deallocate_a(out_fft)
759 case (standard_fourier)
760 call kick_end(kick)
761 safe_deallocate_a(tdrho_a)
762 safe_deallocate_a(wdrho_a)
763 end select
764
765 pop_sub(convert_transform)
766 end subroutine convert_transform
767 ! ---------------------------------------------------------
770 subroutine convert_operate(mesh, namespace, space, ions, mc, outp)
771 class(mesh_t), intent(in) :: mesh
772 type(namespace_t), intent(in) :: namespace
773 class(space_t), intent(in) :: space
774 type(ions_t), intent(in) :: ions
775 type(multicomm_t), intent(in) :: mc
776 type(output_t) , intent(in) :: outp
777
778 integer :: ierr, ip, i_op, length, n_operations, output_i
779 type(block_t) :: blk
780 type(restart_t) :: restart
781 type(unit_t) :: units
782 real(real64) :: f_re, f_im
783 real(real64), allocatable :: tmp_ff(:), scalar_ff(:)
784
785 character(len=200) :: var, scalar_expression
786 character(len=MAX_PATH_LEN) :: folder, filename, out_folder, out_filename
787
788 push_sub(convert_operate)
789
790 !%Variable ConvertScalarOperation
791 !%Type block
792 !%Section Utilities::oct-convert
793 !%Description
794 !% This variable is used to generate a new mesh function as a linear combination
795 !% different mesh function having the same mesh. Each row defines an operation for
796 !% for a single mesh function.
797 !% The format of the block is the following: <br>
798 !% 'variable name' | 'folder' | 'file' | 'operation'
799 !%End
800 ! First, find out if there is a ConvertScalarOperation block.
801 n_operations = 0
802 if (parse_block(namespace, 'ConvertScalarOperation', blk) == 0) then
803 n_operations = parse_block_n(blk)
804 end if
805
806 if (n_operations == 0) then
807 write(message(1),'(a)')'No operations found. Check the input file'
808 call messages_fatal(1)
809 end if
810
811 !%Variable ConvertOutputFolder
812 !%Type string
813 !%Section Utilities::oct-convert
814 !%Description
815 !% The folder name where the output files will be write. The default is
816 !% <tt>convert</tt>.
817 !%End
818 call parse_variable(namespace, 'ConvertOutputFolder', "convert", out_folder)
819 call add_last_slash(out_folder)
820 call io_mkdir(out_folder, namespace)
821
822 !%Variable ConvertOutputFilename
823 !%Type string
824 !%Default "density"
825 !%Section Utilities::oct-convert
826 !%Description
827 !% Output filename. The name of the file in which the converted mesh function will be
828 !% written in the format specified in <tt>OutputFormat</tt>.
829 !%End
830 call parse_variable(namespace, 'ConvertOutputFilename', 'density', out_filename)
831
832 safe_allocate(tmp_ff(1:mesh%np))
833 safe_allocate(scalar_ff(1:mesh%np))
834 scalar_ff = m_zero
835
836 do i_op = 1, n_operations
837 !read variable name
838 call parse_block_string(blk, i_op-1, 0, var)
839 !read folder path
840 call parse_block_string(blk, i_op-1, 1, folder)
841 call add_last_slash(folder)
842 !read file
843 call parse_block_string(blk, i_op-1, 2, filename)
844 ! Delete the extension if present
845 length = len_trim(filename)
846 if (length > 4) then
847 if (filename(length-3:length) == '.obf') then
848 filename = trim(filename(1:length-4))
849 end if
850 end if
851 ! FIXME: why only real functions? Please generalize.
852 ! TODO: check if mesh function are real or complex.
853 call restart%init(namespace, restart_undefined, restart_type_load, mc, ierr, &
854 dir=trim(folder), mesh = mesh, exact=.true.)
855 if (ierr == 0) then
856 call restart%read_mesh_function(space, trim(filename), mesh, tmp_ff, ierr)
857 else
858 write(message(1),'(2a)') "Failed to read from file ", trim(filename)
859 write(message(2), '(2a)') "from folder ", trim(folder)
860 call messages_fatal(2)
861 end if
862 !read scalar expression
863 call parse_block_string(blk, i_op-1, 3, scalar_expression)
864
865 do ip = 1, mesh%np
866 call parse_expression(f_re, f_im, trim(var), real(tmp_ff(ip), real64), trim(scalar_expression))
867 !TODO: implement use of complex functions.
868 scalar_ff(ip) = scalar_ff(ip) + f_re
869 end do
870
871 call restart%end()
872
873 end do
874
875 call parse_block_end(blk)
876
877 call mesh%mpi_grp%barrier()
878 ! Write the corresponding output
879 !TODO: add variable ConvertFunctionType to select the type(density, wfs, potential, ...)
880 ! and units of the conversion.
881 units = units_out%length**(-space%dim)
882 do output_i = lbound(outp%how, 1), ubound(outp%how, 1)
883 if (outp%how(output_i) /= 0) then
884 call dio_function_output(outp%how(output_i), trim(out_folder), trim(out_filename), namespace, space, mesh, &
885 scalar_ff, units, ierr, pos=ions%pos, atoms=ions%atom)
886 end if
887 end do
888
889 safe_deallocate_a(tmp_ff)
890 safe_deallocate_a(scalar_ff)
891
892 pop_sub(convert_operate)
893 end subroutine convert_operate
894end program
895
896!! Local Variables:
897!! mode: f90
898!! coding: utf-8
899!! End:
program oct_convert
This utility runs in parallel and can be used for post-processing of the results of Output.
Definition: convert.F90:118
subroutine convert_low(mesh, namespace, space, ions, psolver, mc, basename, in_folder, c_start, c_end, c_step, outp, iterate_folder, subtract_file, ref_name, ref_folder)
Giving a range of input files, it writes the corresponding output files.
Definition: convert.F90:372
subroutine convert_operate(mesh, namespace, space, ions, mc, outp)
Given a set of mesh function operations it computes a a resulting mesh function from linear combinati...
Definition: convert.F90:866
subroutine convert_transform(mesh, namespace, space, ions, mc, kpoints, basename, in_folder, c_start, c_end, c_step, outp, subtract_file, ref_name, ref_folder)
Giving a range of input files, it computes the Fourier transform of the file.
Definition: convert.F90:503
initialize a batch with existing memory
Definition: batch.F90:277
Each program/utility that needs to use the getopt features should have an interface here – the defini...
Prints out to iunit a message in the form: ["InputVariable" = value] where "InputVariable" is given b...
Definition: messages.F90:182
static void convert(multi *in, multi *out, int t_in, int t_out)
Definition: io_binary.c:5135
This module implements batches of mesh functions.
Definition: batch.F90:135
This module handles the calculation mode.
type(calc_mode_par_t), public calc_mode_par
Singleton instance of parallel calculation mode.
integer, parameter, public p_strategy_states
parallelization in states
subroutine, public getopt_init(ierr)
Initializes the getopt machinery. Must be called before attempting to parse the options....
subroutine, public getopt_end
Fast Fourier Transform module. This module provides a single interface that works with different FFT ...
Definition: fft.F90:120
subroutine, public fft_init(this, nn, dim, type, library, optimize, optimize_parity, comm, mpi_grp, use_aligned)
Definition: fft.F90:411
subroutine, public fft_all_init(namespace)
initialize the table
Definition: fft.F90:280
subroutine, public fft_all_end()
delete all plans
Definition: fft.F90:391
integer, parameter, public fft_real
Definition: fft.F90:180
integer, parameter, public fftlib_fftw
Definition: fft.F90:185
real(real64), parameter, public m_two
Definition: global.F90:192
subroutine, public global_end()
Finalise parser varinfo file, and MPI.
Definition: global.F90:417
real(real64), parameter, public m_zero
Definition: global.F90:190
real(real64), parameter, public m_pi
some mathematical constants
Definition: global.F90:188
subroutine, public global_init(communicator)
Initialise Octopus.
Definition: global.F90:360
subroutine, public dwrite_header(fname, np_global, ierr)
subroutine, public dio_function_output(how, dir, fname, namespace, space, mesh, ff, unit, ierr, pos, atoms, grp, root)
Top-level IO routine for functions defined on the mesh.
Definition: io.F90:116
subroutine, public io_init(defaults)
If the argument defaults is present and set to true, then the routine will not try to read anything f...
Definition: io.F90:165
subroutine, public io_close(iunit, grp)
Definition: io.F90:466
subroutine, public io_end()
Definition: io.F90:270
subroutine, public io_mkdir(fname, namespace, parents)
Definition: io.F90:360
integer function, public io_open(file, namespace, action, status, form, position, die, recl, grp)
Definition: io.F90:401
subroutine, public kick_end(kick)
Definition: kick.F90:796
subroutine, public kick_init(kick, namespace, space, kpoints, nspin)
Definition: kick.F90:225
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:120
subroutine, public messages_end()
Definition: messages.F90:279
subroutine, public messages_print_with_emphasis(msg, iunit, namespace)
Definition: messages.F90:904
character(len=512), private msg
Definition: messages.F90:167
subroutine, public messages_init(output_dir)
Definition: messages.F90:226
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:531
subroutine, public print_date(str)
Definition: messages.F90:989
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:416
subroutine, public messages_experimental(name, namespace)
Definition: messages.F90:1069
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:600
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
type(namespace_t), public global_namespace
Definition: namespace.F90:134
this module contains the low-level part of the output system
Definition: output_low.F90:117
this module contains the output system
Definition: output.F90:117
subroutine, public parser_init()
Initialise the Octopus parser.
Definition: parser.F90:453
subroutine, public parser_end()
End the Octopus parser.
Definition: parser.F90:484
integer function, public parse_block(namespace, name, blk, check_varinfo_)
Definition: parser.F90:621
subroutine, public dpoisson_solve(this, namespace, pot, rho, all_nodes, kernel, reset)
Calculates the Poisson equation. Given the density returns the corresponding potential.
Definition: poisson.F90:871
subroutine, public profiling_end(namespace)
Definition: profiling.F90:415
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 profiling_init(namespace)
Create profiling subdirectory.
Definition: profiling.F90:257
integer, parameter, public restart_undefined
Definition: restart.F90:156
integer, parameter, public restart_type_load
Definition: restart.F90:183
subroutine, public spectrum_fourier_transform(method, transform, noise, time_start, time_end, t0, time_step, time_function, energy_start, energy_end, energy_step, energy_function)
Computes the sine, cosine, (or "exponential") Fourier transform of the real function given in the tim...
Definition: spectrum.F90:2639
subroutine, public spectrum_init(spectrum, namespace, default_energy_step, default_max_energy)
Definition: spectrum.F90:215
integer default
Definition: spectrum.F90:209
subroutine, public spectrum_signal_damp(damp_type, damp_factor, time_start, time_end, t0, time_step, time_function)
Definition: spectrum.F90:2553
subroutine, public add_last_slash(str)
Adds a '/' in the end of the string, only if it missing. Useful for directories.
Definition: string.F90:163
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
Definition: unit.F90:134
character(len=20) pure function, public units_abbrev(this)
Definition: unit.F90:225
This module defines the unit system, used for input and output.
type(unit_system_t), public units_out
subroutine, public unit_system_init(namespace)
type(unit_system_t), public units_inp
the units systems for reading and writing
This module is intended to contain simple general-purpose utility functions and procedures.
Definition: utils.F90:120
character(len=256) function, public get_config_opts()
Definition: utils.F90:365
character(len=256) function, public get_optional_libraries()
Definition: utils.F90:371
subroutine, public print_header()
This subroutine prints the logo followed by information about the compilation and the system....
Definition: utils.F90:304
Class defining batches of mesh functions.
Definition: batch.F90:161
Class describing the electron system.
Definition: electrons.F90:220
Describes mesh distribution to nodes.
Definition: mesh.F90:187
Stores all communicators and groups.
Definition: multicomm.F90:208
output handler class
Definition: output_low.F90:166
int true(void)