io_function.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module io_function_oct_m
  use affine_coordinates_oct_m
  use atom_oct_m
  use box_oct_m
  use box_parallelepiped_oct_m
  use comm_oct_m
  use cube_function_oct_m
  use cube_oct_m
  use debug_oct_m
  use distributed_oct_m
  use fft_oct_m
  use global_oct_m
  use index_oct_m
  use io_oct_m
  use io_binary_oct_m
  use io_csv_oct_m
  use iso_c_binding
  use kpoints_oct_m
  use lattice_vectors_oct_m
  use mesh_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use mpi_oct_m
  use mpi_debug_oct_m
  use namespace_oct_m
#if defined(HAVE_NETCDF)
  use netcdf
#endif
  use par_vec_oct_m
  use parser_oct_m
  use partition_oct_m
  use profiling_oct_m
  use space_oct_m
  use species_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use utils_oct_m
  use varinfo_oct_m
  use vtk_oct_m
 
  implicit none
 
  private
  public ::                         &
    io_function_read_what_how_when, &
    io_function_fill_how,           &
    write_canonicalized_xyz_file,   &
    write_standard_xyz_file,        &
    write_xsf_geometry,             &
    write_xsf_geometry_file,        &
    dio_function_input,             &
    zio_function_input,             &
    dio_function_output,            &
    zio_function_output,            &
    io_function_output_vector,      &
    dio_function_output_global,     &
    zio_function_output_global,     &
    io_function_output_supercell,   &
    io_function_output_vector_bz,   &
    io_function_output_global_bz
 
#if defined(HAVE_NETCDF)
  public ::                         &
    dout_cf_netcdf,                 &
    zout_cf_netcdf
#endif
 
  integer, parameter, private ::    &
    DOUTPUT_KIND      =    1,       &
    zoutput_kind      =   -1
 
  character(len=3), parameter ::    &
    index2label(3) = (/ 're ', 'im ', 'abs' /)
 
 
  interface io_function_output_vector
    module procedure dio_function_output_vector, zio_function_output_vector
  end interface io_function_output_vector
 
  interface io_function_output_vector_bz
    module procedure dio_function_output_vector_bz, zio_function_output_vector_bz
  end interface io_function_output_vector_bz
 
  interface io_function_output_global_bz
    module procedure dio_function_output_global_bz, zio_function_output_global_bz
  end interface io_function_output_global_bz
 
  interface io_function_output_supercell
    module procedure dio_function_output_supercell, zio_function_output_supercell
  end interface io_function_output_supercell
 
 
contains
 
  ! ---------------------------------------------------------
  subroutine io_function_read_what_how_when(namespace, space, what, how, output_interval, &
    what_tag_in, how_tag_in, output_interval_tag_in, ignore_error)
    type(namespace_t), intent(in)           :: namespace
    class(space_t),    intent(in)           :: space
    logical,           intent(inout)        :: what(MAX_OUTPUT_TYPES)
    integer(int64),       intent(out)          :: how(0:MAX_OUTPUT_TYPES)
    integer,           intent(out)          :: output_interval(0:MAX_OUTPUT_TYPES)
    character(len=*),  optional, intent(in) :: what_tag_in
    character(len=*),  optional, intent(in) :: how_tag_in
    character(len=*),  optional, intent(in) :: output_interval_tag_in
    logical, optional, intent(in)           :: ignore_error
 
    type(block_t) :: blk
    integer :: ncols, nrows, iout, column_index, what_i
    integer(int64) :: what_no_how(13)
    character(len=80) :: what_tag
    character(len=80) :: how_tag
    character(len=80) :: output_interval_tag
    character(len=80) :: output_column_marker
    push_sub(io_function_read_what_how_when)
 
    how = 0_8
    output_interval = 0_8
    ncols = 0
    what_tag = optional_default(what_tag_in, 'Output')
    how_tag = optional_default(how_tag_in, 'OutputFormat')
    output_interval_tag = optional_default(output_interval_tag_in, 'OutputInterval')
 
    call messages_obsolete_variable(namespace, 'OutputHow', 'OutputFormat')
    call messages_obsolete_variable(namespace, 'Output_KPT', 'Output')
    call messages_obsolete_variable(namespace, 'OutputLDA_U' , 'Output')
    call messages_obsolete_variable(namespace, 'OutputEvery', 'OutputInterval')
 
    !%Variable Output
    !%Type block
    !%Default none
    !%Section Output
    !%Description
    !% Specifies what to print.
    !% Each output must be in a separate row. Optionally individual output formats and output intervals can be defined
    !% for each row or they can be read separately from <tt>OutputFormat</tt> and <tt>OutputInterval</tt> variables
    !% in the input file.
    !% The output files are written at the end of the run into the output directory for the
    !% relevant kind of run (<i>e.g.</i> <tt>static</tt> for <tt>CalculationMode = gs</tt>).
    !% Time-dependent simulations print only per iteration, including always the last. The frequency of output per iteration
    !% (available for <tt>CalculationMode</tt> = <tt>gs</tt>, <tt>unocc</tt>,  <tt>td</tt>, and <tt>opt_control</tt>)
    !% is set by <tt>OutputInterval</tt> and the directory is set by <tt>OutputIterDir</tt>.
    !% For linear-response run modes, the derivatives of many quantities can be printed, as listed in
    !% the options below. Indices in the filename are labelled as follows:
    !% <tt>sp</tt> = spin (or spinor component), <tt>k</tt> = <i>k</i>-point, <tt>st</tt> = state/band.
    !% There is no tag for directions, given as a letter. The perturbation direction is always
    !% the last direction for linear-response quantities, and a following +/- indicates the sign of the frequency.
    !%
    !% Example (minimal):
    !% <br><br><tt>%Output
    !% <br>&nbsp;&nbsp;density
    !% <br>&nbsp;&nbsp;potential
    !% <br>%<br></tt>
    !%
    !% Example (with OutputFormat):
    !% <br><br><tt>%Output
    !% <br>&nbsp;&nbsp;density   | cube + axis_z
    !% <br>&nbsp;&nbsp;potential | cube
    !% <br>%<br></tt>
    !%
    !% Example (with OutputFormat, incomplete):
    !% <br><br><tt>%Output
    !% <br>&nbsp;&nbsp;density   | cube + axis_z
    !% <br>&nbsp;&nbsp;potential
    !% <br>%<br></tt>
    !%
    !% Example (tagged):
    !% <br><br><tt>%Output
    !% <br>&nbsp;&nbsp;density   | "output_format" | cube + axis_z | "output_interval" | 50
    !% <br>&nbsp;&nbsp;potential | "output_format" | cube          | "output_interval" | 20
    !% <br>%<br></tt>
    !%
    !% Example (tagged, incomplete):
    !% <br><br><tt>%Output
    !% <br>&nbsp;&nbsp;density   | "output_format"   | cube + axis_z
    !% <br>&nbsp;&nbsp;potential | "output_interval" | 20
    !% <br>%<br></tt>
    !% Missing information for the incomplete blocks will be parsed form the out-of-block
    !% definitions. It is also possible to mix the order of columns in the tagged format.
    !% See <tt>OutputFormat</tt>, and <tt>OutputInterval</tt>.
    !%Option potential 1
    !% Outputs Kohn-Sham potential, separated by parts. File names are <tt>v0</tt> for
    !% the local part of the ionic potential, <tt>vc</tt> for the classical potential (if it exists),
    !% <tt>vh</tt> for the Hartree potential, <tt>vks</tt> for the local part of the Kohn-Sham potential, and
    !% <tt>vxc-</tt> for the exchange-correlation potentials. For <tt>vks</tt> and <tt>vxc</tt>,
    !% a suffix for spin is added in the spin-polarized case.
    !%Option density 2
    !% Outputs density. The output file is called <tt>density-</tt>, or <tt>lr_density-</tt> in linear response.
    !%Option wfs 3
    !% Outputs wavefunctions. Which wavefunctions are to be printed is specified
    !% by the variable <tt>OutputWfsNumber</tt> -- see below. The output file is called
    !% <tt>wf-</tt>, or <tt>lr_wf-</tt> in linear response.
    !%Option wfs_sqmod 4
    !% Outputs modulus squared of the wavefunctions.
    !% The output file is called <tt>sqm-wf-</tt>. For linear response, the filename is <tt>sqm_lr_wf-</tt>.
    !%Option geometry 5
    !% Outputs file containing the coordinates of the atoms treated within quantum mechanics.
    !% If <tt>OutputFormat = xyz</tt>, the file is called <tt>geometry.xyz</tt>; a
    !% file <tt>crystal.xyz</tt> is written with a supercell geometry if the system is periodic;
    !% if point charges were defined in the PDB file (see <tt>PDBCoordinates</tt>), they will be output
    !% in the file <tt>geometry_classical.xyz</tt>.
    !% If <tt>OutputFormat = xcrysden</tt>, a file called <tt>geometry.xsf</tt> is written.
    !%Option current 6
    !% Outputs the total current density. The output file is called <tt>current-</tt>.
    !% For linear response, the filename is <tt>lr_current-</tt>.
    !%Option ELF 7
    !% Outputs electron localization function (ELF). The output file is called <tt>elf-</tt>.
    !% Only in 2D and 3D.
    !%Option ELF_basins 8
    !% Outputs basins of attraction of the ELF. The output file is called
    !% <tt>elf_rs_basins.info</tt>. Only in 2D and 3D.
    !%Option Bader 9
    !% Outputs Laplacian of the density which shows lone pairs, bonded charge concentrations
    !% and regions subject to electrophilic or nucleophilic attack.
    !% See RF Bader, <i>Atoms in Molecules: A Quantum Theory</i> (Oxford Univ. Press, Oxford, 1990).
    !%Option el_pressure 10
    !% Outputs electronic pressure. See Tao, Vignale, and Tokatly, <i>Phys Rev Lett</i> <b>100</b>, 206405 (2008).
    !%Option matrix_elements 11
    !% Outputs a series of matrix elements of the Kohn-Sham states. What is output can
    !% be controlled by the <tt>OutputMatrixElements</tt> variable.
    !%Option pol_density 12
    !% Outputs dipole-moment density <tt>dipole_density-</tt>, or polarizability density <tt>alpha_density-</tt>
    !% in linear response. If <tt>ResponseMethod = finite_differences</tt>, the hyperpolarizability density
    !% <tt>beta_density-</tt> is also printed.
    !%Option mesh_r 13
    !% Outputs values of the coordinates over the grid. Files
    !% will be called <tt>mesh_r-</tt> followed by the direction.
    !%Option kinetic_energy_density 14
    !% Outputs kinetic-energy density, defined as:
    !%
    !% <math>\tau_\sigma(\vec{r}) = \sum_{i=1}^{N_\sigma}
    !%  \left| \vec{\nabla} \phi_{i\sigma}(\vec{r}) \right|^2\,. </math>
    !%
    !% The index <math>\sigma</math> is the spin index for the spin-polarized case,
    !% or if you are using spinors. For spin-unpolarized calculations, you
    !% get the total kinetic-energy density. The previous expression assumes full
    !% or null occupations. If fractional occupation numbers, each term in the sum
    !% is weighted by the occupation. Also, if we are working with an infinite
    !% system, all <i>k</i>-points are summed up, with their corresponding weights. The
    !% files will be called <tt>tau-sp1</tt> and <tt>tau-sp2</tt>, if the spin-resolved kinetic
    !% energy density is produced (runs in spin-polarized and spinors mode), or
    !% only <tt>tau</tt> if the run is in spin-unpolarized mode.
    !%
    !% Note that the code does not the U(1)gauge-invariant kinetic energy density, but simply the kinetic energy density.
    !%Option dos 15
    !% Outputs density of states. See <tt>DOSEnergyMax</tt>, <tt>DOSEnergyMin</tt>, <tt>DOSEnergyPoints</tt>,
    !% and <tt>DOSGamma</tt>.
    !%Option tpa 16
    !% Outputs transition-potential approximation (TPA) matrix elements, using <math>\vec{q}</math>-vector specified
    !% by <tt>MomentumTransfer</tt>.
    !%Option forces 17
    !% Outputs file <tt>forces.xsf</tt> containing structure and forces on the atoms as
    !% a vector associated with each atom, which can be visualized with XCrySDen.
    !%Option wfs_fourier 18
    !% (Experimental) Outputs wavefunctions in Fourier space. This is
    !% only implemented for the ETSF file format output. The file will
    !% be called <tt>wfs-pw-etsf.nc</tt>.
    !%Option xc_density 19
    !% Outputs the XC density, which is the charge density that
    !% generates the XC potential. (This is <math>-1/4\pi</math> times
    !% the Laplacian of the XC potential). The files are called <tt>nxc</tt>.
    !%Option PES_wfs 20
    !% Outputs the photoelectron wavefunctions. The file name is <tt>pes_wfs-</tt>
    !% plus the orbital number.
    !%Option PES_density 21
    !% Outputs the photolectron density. Output file is <tt>pes_dens-</tt> plus spin species if
    !% spin-polarized calculation is performed.
    !%Option PES 22
    !% Outputs the time-dependent photoelectron spectrum.
    !%Option BerkeleyGW 23
    !% Output for a run with <a href=http://www.berkeleygw.org>BerkeleyGW</a>.
    !% See <tt>Output::BerkeleyGW</tt> for further specification.
    !%Option delta_perturbation 24
    !% Outputs the "kick", or time-delta perturbation applied to compute optical response in real time.
    !%Option external_td_potential 25
    !% Outputs the (scalar) time-dependent potential.
    !%Option mmb_wfs 26
    !% Triggers the ModelMB wavefunctions to be output for each state.
    !%Option mmb_den 27
    !% Triggers the ModelMB density matrix to be output for each state, and the particles
    !% specified by the <tt>DensitytoCalc</tt> block. Calculates, and outputs, the reduced density
    !% matrix. For the moment the trace is made over the second dimension, and
    !% the code is limited to 2D. The idea is to model <i>N</i> particles in 1D as an
    !% <i>N</i>-dimensional non-interacting problem, then to trace out <i>N</i>-1 coordinates.
    !%Option potential_gradient 28
    !% Prints the gradient of the potential.
    !%Option energy_density 29
    !% Outputs the total energy density to a file called
    !% <tt>energy_density</tt>.
    !%Option heat_current 30
    !% Outputs the total heat current density. The output file is
    !% called <tt>heat_current-</tt>.
    !%Option photon_correlator 31
    !% Outputs the electron-photon correlation function. The output file is
    !% called <tt>photon_correlator</tt>.
    !%Option J_flow 32
    !% todo: document J_flow option!
    !%Option current_kpt 33
    !% Outputs the current density resolved in momentum space. The output file is called <tt>current_kpt-</tt>.
    !%Option density_kpt 34
    !% Outputs the electronic density resolved in momentum space.
    !%Option occ_matrices 35
    !% Only for DFT+U calculations.
    !% Outputs the occupation matrices of DFT+U
    !%Option effectiveU 36
    !% Only for DFT+U calculations.
    !% Outputs the value of the effectiveU for each atoms
    !%Option magnetization 37
    !% Only for DFT+U calculations.
    !% Outputs file containing structure and magnetization of the localized subspace
    !% on the atoms as a vector associated with each atom, which can be visualized.
    !% For the moment, it only works if a +U is added on one type of orbital per atom.
    !%Option local_orbitals 38
    !% Only for DFT+U calculations.
    !% Outputs the localized orbitals that form the correlated subspace
    !%Option kanamoriU 39
    !% Only for DFT+U calculations.
    !% Outputs the Kanamori interaction parameters U, U`, and J.
    !% These parameters are not determined self-consistently, but are taken from the
    !% occupation matrices and Coulomb integrals comming from a standard +U calculation.
    !%Option xc_torque 40
    !% Outputs the exchange-correlation torque. Only for the spinor case and in the 3D case.
    !%Option eigenval_kpt 41
    !% Outputs the eigenvalues resolved in momentum space, with one file for each band.
    !%Option stress 42
    !% Outputs the stress tensor and each of its contributing terms
    !%Option current_dia 43
    !% Outputs the diamagnetic current density from a non-uniform vector potential.
    !% The output file is called <tt>current_dia-</tt>.
    !%Option jdos 44
    !% Outputs the joint density of states.
    !% The same variables as for the regular DOS are used to control the energies and broadening.
    !%Option ldos 45
    !% Outputs the local density of states
    !% The broadening uses the same value as the regular DOS. The energies at which the LDOS is computed are
    !% specified by the <tt>LDOSEnergies</tt> block.
    !%Option pdos 46
    !% Outputs the projected density of states.
    !% At the moment, the PDOS is computed from the bare pseudo-atomic orbitals, directly taken from
    !% the pseudopotentials. The orbitals are not orthonormalized, in order to preserve their
    !% atomic orbitals character. As a consequence, the sum of the different PDOS does not integrate
    !% to the total DOS.
    !%
    !% The radii of the orbitals are controled by the threshold defined by <tt>AOThreshold</tt>,
    !% and the fact that they are normalized or not by <tt>AONormalize</tt>.
    !%End
 
    !%Variable OutputFormat
    !%Type flag
    !%Default 0
    !%Section Output
    !%Description
    !% Describes the format of the output files.
    !% This variable can also be defined inside the <tt>Output</tt> block.
    !% See <tt>Output</tt>.
    !% Example: <tt>axis_x + plane_x + dx</tt>
    !%Option axis_x bit(0)
    !% The values of the function on the <i>x</i> axis are printed. The string <tt>.y=0,z=0</tt> is appended
    !% to previous file names.
    !% For non-orthogonal cells, this correspond to the first lattice vector instead the x axis.
    !%Option axis_y bit(1)
    !% The values of the function on the <i>y</i> axis are printed. The string <tt>.x=0,z=0</tt> is appended
    !% to previous file names.
    !% For non-orthogonal cells, this correspond to the second lattice vector instead the y axis.
    !%Option axis_z bit(2)
    !% The values of the function on the <i>z</i> axis are printed. The string <tt>.x=0,y=0</tt> is appended
    !% to previous file names.
    !% For non-orthogonal cells, this correspond to the third lattice vector instead the z axis.
    !%Option plane_x bit(3)
    !% A plane slice at <i>x</i> = 0 is printed. The string <tt>.x=0</tt> is appended
    !% to previous file names.
    !%Option plane_y bit(4)
    !% A plane slice at <i>y</i> = 0 is printed. The string <tt>.y=0</tt> is appended
    !% to previous file names.
    !%Option plane_z bit(5)
    !% A plane slice at <i>z</i> = 0 is printed. The string <tt>.z=0</tt> is appended to
    !% previous file names.
    !%Option dx bit(6)
    !% For printing three-dimensional information, the open-source program
    !% visualization tool <a href=http://www.opendx.org>OpenDX</a> can be used. The string
    !% <tt>.dx</tt> is appended to previous file names. Available only in 3D.
    !%Option netcdf bit(7)
    !% Outputs in <a href=http://www.unidata.ucar.edu/packages/netcdf>NetCDF</a> format. This file
    !% can then be read, for example, by OpenDX. The string <tt>.ncdf</tt> is appended to previous file names.
    !% Requires the NetCDF library. Only writes the real part of complex functions.
    !%Option mesh_index bit(8)
    !% Generates output files of a given quantity (density, wavefunctions, ...) which include
    !% the internal numbering of mesh points. Since this mode produces large datafiles this is only
    !% useful for small meshes and debugging purposes.
    !% The output can also be used to display the mesh directly. A Gnuplot script for mesh visualization
    !% can be found under <tt>PREFIX/share/octopus/util/display_mesh_index.gp</tt>.
    !%Option xcrysden bit(9)
    !% A format for printing structures and three-dimensional information, which can be visualized by
    !% the free open-source program <a href=http://www.xcrysden.org>XCrySDen</a> and others. The string
    !% <tt>.xsf</tt> is appended to previous file names. Note that lattice vectors and coordinates are as
    !% specified by <tt>UnitsOutput</tt>. Available in 2D and 3D.
    !%Option matlab bit(10)
    !% In combination with <tt>plane_x</tt>, <tt>plane_y</tt> and
    !% <tt>plane_z</tt>, this option produces output files which are
    !% suitable for 2D Matlab functions like <tt>mesh()</tt>,
    !% <tt>surf()</tt>, or <tt>waterfall()</tt>. To load these files
    !% into Matlab you can use, <i>e.g.</i>
    !%<tt>
    !%   >> density = load('static/density-1.x=0.matlab.abs');
    !%   >> mesh(density);
    !%</tt>
    !%Option meshgrid bit(11)
    !% Outputs in Matlab mode the internal mesh in a format similar to
    !%<tt>
    !%   >> [x,y] = meshgrid(-2:.2:2,-1:.15:1)
    !%</tt>
    !% The <i>x</i> meshgrid is contained in a file <tt>*.meshgrid.x</tt> and the <i>y</i>-grid can be found in
    !% <tt>*.meshgrid.y</tt>.
    !%Option boundary_points bit(12)
    !% This option includes the output of the mesh enlargement. Default is without.
    !% Supported only by <tt>binary</tt>, <tt>axis</tt>, <tt>plane</tt>, <tt>mesh_index</tt>,
    !% and <tt>matlab</tt> formats.
    !% Not all types of <tt>Output</tt> will have this information available. Not supported when parallel in domains.
    !%Option binary bit(13)
    !% Plain binary, new format.
    !%Option xyz bit(15)
    !% Geometry will be output in XYZ format. Does not affect other outputs.
    !%Option cube bit(16)
    !% Generates output in the <a href=http://paulbourke.net/dataformats/cube>cube file format</a>.
    !% Available only in 3D. Only writes the real part of complex functions.
    !% This output format always uses atomic units.
    !%Option bild bit(19)
    !% Generates output in <a href=http://plato.cgl.ucsf.edu/chimera/docs/UsersGuide/bild.html>BILD format</a>.
    !%Option vtk bit(20)
    !% Generates output in <a href=http://www.vtk.org/VTK/img/file-formats.pdf>VTK legacy format</a>.
    !%Option integrate_xy bit(21)
    !% Integrates the function in the x-y plane and the result on the <i>z</i> axis is printed.
    !%Option integrate_xz bit(22)
    !% Integrates the function in the x-z plane and the result on the <i>y</i> axis is printed
    !%Option integrate_yz bit(23)
    !% Integrates the function in the y-z plane and the result on the <i>x</i> axis is printed
    !%Option ascii bit(24)
    !% Plain text format regardless of dimensionality. For the moment only employed by the oct-phototoelectron_spectrum
    !% post-processing utility.
    !%Option poscar bit(25)
    !% Write geometry in POSCAR format.
    !% Note that atoms are not ordered by labels and that species names are taken from the Octopus input,
    !% and are not truncated to two letters.
    !%End
 
    !%Variable OutputInterval
    !%Type integer
    !%Default 50
    !%Section Output
    !%Description
    !% The output requested by variable <tt>Output</tt> is written
    !% to the directory <tt>OutputIterDir</tt>
    !% when the iteration number is a multiple of the <tt>OutputInterval</tt> variable.
    !% Subdirectories are named Y.X, where Y is <tt>td</tt>, <tt>scf</tt>, or <tt>unocc</tt>, and
    !% X is the iteration number. To use the working directory, specify <tt>"."</tt>
    !% Must be >= 0. If it is 0, then no output is written. For <tt>gs</tt> and <tt>unocc</tt>
    !% calculations, <tt>OutputDuringSCF</tt> must be set too for this output to be produced.
    !% This variable can also be defined inside the <tt>Output</tt> block.
    !% See <tt>Output</tt>.
    !%End
 
 
    what_no_how = (/ &
      option__output__matrix_elements, &
      option__output__berkeleygw, &
      option__output__dos, &
      option__output__tpa, &
      option__output__mmb_den, &
      option__output__j_flow, &
      option__output__occ_matrices, &
      option__output__effectiveu, &
      option__output__magnetization, &
      option__output__kanamoriu, &
      option__output__stress, &
      option__output__jdos, &
      option__output__pdos &
      /)
 
    if (parse_block(namespace, what_tag, blk) == 0) then
      nrows = parse_block_n(blk)
      do iout = 0, nrows - 1
        ncols = max(ncols , parse_block_cols(blk, iout))
      end do
 
      if (ncols == 1) then
        !new format, Type 0
        !%Output
        !  density
        !  wfs
        !%
        do iout = 1, nrows
          call parse_block_integer(blk, iout - 1, 0, what_i)
          if (.not. varinfo_valid_option(what_tag, what_i)) then
            call messages_input_error(namespace, what_tag)
          end if
          if (what_i > 0) then
            what(what_i) = .true.
            call parse_variable(namespace, output_interval_tag, 50, output_interval(what_i))
            if (((what_tag == 'Output') .and. (.not. any(what_no_how == what_i)))&
              .or. (what_tag /= 'Output')) then
              call parse_variable(namespace, how_tag, 0, how(what_i))
              if (.not. varinfo_valid_option(how_tag, how(what_i), is_flag=.true.)) then
                call messages_input_error(namespace, how_tag)
              end if
            end if
          end if
        end do
      else if (ncols == 2) then
        !new format, Type 1
        !%Output
        !  density | cube + axis_z
        !  wfs     | cube
        !%
 
        do iout = 1, nrows
          call parse_block_integer(blk, iout - 1, 0, what_i)
          if (.not. varinfo_valid_option(what_tag, what_i)) then
            call messages_input_error(namespace, what_tag)
          end if
          if (what_i > 0) then
            what(what_i) = .true.
            call parse_variable(namespace, output_interval_tag, 50, output_interval(what_i))
            if (((what_tag == 'Output') .and. (.not. any(what_no_how == what_i)))&
              .or. (what_tag /= 'Output')) then
              call parse_block_integer(blk, iout - 1, 1, how(what_i))
              if (how(what_i) == 0) call parse_variable(namespace, how_tag, 0, how(what_i))
              if (.not. varinfo_valid_option(how_tag, how(what_i), is_flag=.true.)) then
                call messages_input_error(namespace, how_tag)
              end if
            end if
          end if
        end do
 
      else
        !new format, Type 2 (tagged)
        !%Output
        !  density | "output_interval" | 10   | "output_format"   | cube + axis_z
        !  wfs     | "output_format"   | cube | "output_interval" | 50
        !%
        !
        ! OR
        !
        !%Output
        !  density | "output_interval" | 10   | "output_format" | cube + axis_z
        !  wfs     | "output_format"   | cube
        !%
        do iout = 1, nrows
          call parse_block_integer(blk, iout - 1, 0, what_i)
          if (.not. varinfo_valid_option(what_tag, what_i)) then
            call messages_input_error(namespace, what_tag)
          end if
          if (what_i > 0) then
            what(what_i) = .true.
            do column_index = 0, 1
              call parse_block_string(blk, iout - 1, 1 + column_index * 2, output_column_marker)
              if (output_column_marker == 'output_interval') then
                call parse_block_integer(blk, iout - 1, 2 + column_index * 2, output_interval(what_i))
              else if (output_column_marker == 'output_format') then
                if (((what_tag == 'Output') .and. (.not. any(what_no_how == what_i)))&
                  .or. (what_tag /= 'Output')) then
                  call parse_block_integer(blk, iout - 1, 2 + column_index * 2, how(what_i))
                  if (.not. varinfo_valid_option(how_tag, how(what_i), is_flag=.true.)) then
                    call messages_input_error(namespace, how_tag)
                  end if
                end if
              else if (len_trim(output_column_marker) /= 0) then
                ! Unknown output_column_marker
                call messages_input_error(namespace, what_tag)
              else
                ! no output_column_marker -> full output info is not in this block
                if (output_interval(what_i) == 0) then
                  call parse_variable(namespace, output_interval_tag, 50, output_interval(what_i))
                end if
                if (how(what_i) == 0) then
                  if (((what_tag == 'Output') .and. (.not. any(what_no_how == what_i)))&
                    .or. (what_tag /= 'Output')) then
                    call parse_variable(namespace, how_tag, 0, how(what_i))
                    if (.not. varinfo_valid_option(how_tag, how(what_i), is_flag=.true.)) then
                      call messages_input_error(namespace, how_tag)
                    end if
                  end if
                end if
              end if
            end do
          end if
        end do
      end if
      call parse_block_end(blk)
    else
 
      call messages_variable_is_block(namespace, what_tag)
 
      ! Output block does not exist but we may have OutputHow/OutputInterval
      call parse_variable(namespace, how_tag, 0, how(0))
      call parse_variable(namespace, output_interval_tag, 50, output_interval(0))
    end if
 
 
    do what_i = lbound(what, 1), ubound(what, 1)
      if (what_tag == 'Output') then
        if (what(what_i) .and. (.not. any(what_no_how == what_i))) then
          if (.not. varinfo_valid_option(how_tag, how(what_i), is_flag=.true.)) then
            call messages_input_error(namespace, how_tag)
          end if
 
          if (how(what_i) == 0 .and. .not. optional_default(ignore_error, .false.)) then
            write(message(1), '(a)') 'Must specify output method with variable OutputFormat.'
            call messages_fatal(1, only_root_writes = .true., namespace=namespace)
          end if
 
          ! some modes are not available in some circumstances
          if (space%dim == 1) then
            if (bitand(how(what_i), option__outputformat__axis_y) /= 0) then
              message(1) = "OutputFormat = axis_y not available with Dimensions = 1."
              call messages_fatal(1, namespace=namespace)
            end if
            if (bitand(how(what_i), option__outputformat__plane_z) /= 0) then
              message(1) = "OutputFormat = plane_z not available with Dimensions = 1."
              call messages_fatal(1, namespace=namespace)
            end if
            if (bitand(how(what_i), option__outputformat__xcrysden) /= 0) then
              message(1) = "OutputFormat = xcrysden not available with Dimensions = 1."
              call messages_fatal(1, namespace=namespace)
            end if
          end if
 
          if (space%dim <= 2) then
            if (bitand(how(what_i), option__outputformat__axis_z) /= 0) then
              message(1) = "OutputFormat = axis_z not available with Dimensions <= 2."
              call messages_fatal(1, namespace=namespace)
            end if
            if (bitand(how(what_i), option__outputformat__plane_x) /= 0) then
              message(1) = "OutputFormat = plane_x not available with Dimensions <= 2."
              call messages_fatal(1, namespace=namespace)
            end if
            if (bitand(how(what_i), option__outputformat__plane_y) /= 0) then
              message(1) = "OutputFormat = plane_y not available with Dimensions <= 2."
              call messages_fatal(1, namespace=namespace)
            end if
            if (bitand(how(what_i), option__outputformat__integrate_xy) /= 0) then
              message(1) = "OutputFormat = integrate_xy not available with Dimensions <= 2."
              call messages_fatal(1, namespace=namespace)
            end if
            if (bitand(how(what_i), option__outputformat__integrate_xz) /= 0) then
              message(1) = "OutputFormat = integrate_xz not available with Dimensions <= 2."
              call messages_fatal(1, namespace=namespace)
            end if
            if (bitand(how(what_i), option__outputformat__integrate_yz) /= 0) then
              message(1) = "OutputFormat = integrate_yz not available with Dimensions <= 2."
              call messages_fatal(1, namespace=namespace)
            end if
            if (bitand(how(what_i), option__outputformat__dx) /= 0) then
              message(1) = "OutputFormat = dx not available with Dimensions <= 2."
              call messages_fatal(1, namespace=namespace)
            end if
            if (bitand(how(what_i), option__outputformat__cube) /= 0) then
              message(1) = "OutputFormat = cube not available with Dimensions <= 2."
              call messages_fatal(1, namespace=namespace)
            end if
            if (bitand(how(what_i), option__outputformat__poscar) /= 0) then
              message(1) = "OutputFormat = POSCAR not available with Dimensions <= 2."
              call messages_fatal(1, namespace=namespace)
            end if
          end if
 
          if (space%periodic_dim<=2) then
            if (bitand(how(what_i), option__outputformat__poscar) /= 0) then
              message(1) = "OutputFormat = POSCAR not available with periodic_dim <= 2."
              call messages_fatal(1, namespace=namespace)
            end if
          end if
 
#if !defined(HAVE_NETCDF)
          if (bitand(how(what_i), option__outputformat__netcdf) /= 0) then
            message(1) = 'Octopus was compiled without NetCDF support.'
            message(2) = 'It is not possible to write output in NetCDF format.'
            call messages_fatal(2, namespace=namespace)
          end if
#endif
 
 
        end if
      end if
 
      if (output_interval(what_i) < 0) then
        message(1) = "OutputInterval must be >= 0."
        call messages_fatal(1, namespace=namespace)
      end if
    end do
 
    pop_sub(io_function_read_what_how_when)
  end subroutine io_function_read_what_how_when
 
  ! -------------------------------------------------------------------
  !                       ".", "func", mesh, sb, func, M_ONE, ierr)
  ! -------------------------------------------------------------------
  integer(int64) function io_function_fill_how(where) result(how)
    character(len=*), intent(in) :: where
 
    push_sub(io_function_fill_how)
 
    how = 0
    if (index(where, "AxisX")     /= 0) how = ior(how, option__outputformat__axis_x)
    if (index(where, "AxisY")     /= 0) how = ior(how, option__outputformat__axis_y)
    if (index(where, "AxisZ")     /= 0) how = ior(how, option__outputformat__axis_z)
    if (index(where, "PlaneX")    /= 0) how = ior(how, option__outputformat__plane_x)
    if (index(where, "PlaneY")    /= 0) how = ior(how, option__outputformat__plane_y)
    if (index(where, "PlaneZ")    /= 0) how = ior(how, option__outputformat__plane_z)
    if (index(where, "IntegrateXY") /= 0) how = ior(how, option__outputformat__integrate_xy)
    if (index(where, "IntegrateXZ") /= 0) how = ior(how, option__outputformat__integrate_xz)
    if (index(where, "IntegrateYZ") /= 0) how = ior(how, option__outputformat__integrate_yz)
    if (index(where, "DX")        /= 0) how = ior(how, option__outputformat__dx)
    if (index(where, "XCrySDen")  /= 0) how = ior(how, option__outputformat__xcrysden)
    if (index(where, "Binary")    /= 0) how = ior(how, option__outputformat__binary)
    if (index(where, "MeshIndex") /= 0) how = ior(how, option__outputformat__mesh_index)
    if (index(where, "XYZ")       /= 0) how = ior(how, option__outputformat__xyz)
#if defined(HAVE_NETCDF)
    if (index(where, "NETCDF")    /= 0) how = ior(how, option__outputformat__netcdf)
#endif
    if (index(where, "Cube")      /= 0) how = ior(how, option__outputformat__cube)
    if (index(where, "VTK")       /= 0) how = ior(how, option__outputformat__vtk)
    if (index(where, "ASCII")     /= 0) how = ior(how, option__outputformat__ascii)
    if (index(where, "Bild")      /= 0) how = ior(how, option__outputformat__bild)
 
    pop_sub(io_function_fill_how)
  end function io_function_fill_how
 
  subroutine write_standard_xyz_file(namespace, fname, pos, species, header)
    type(namespace_t),           intent(in) :: namespace
    character(len=*),            intent(in) :: fname
    real(real64),                intent(in) :: pos(:,:)
    character(len=*),            intent(in) :: species(:)
    character(len=*),  optional, intent(in) :: header
 
    integer      :: iunit, n_atoms, iatom, min_dim
    real(real64) :: position(3)
 
    push_sub(write_standard_xyz_file)
 
    n_atoms = size(pos, 2)
    assert(size(species) == n_atoms)
 
    iunit = io_open(trim(fname)//'.xyz', namespace, action='write', position='asis')
 
    write(iunit, '(i6)') n_atoms
    if (present(header)) then
      write(iunit, '(a)') trim(adjustl(header))
    else
      write(iunit, '(a)')
    endif
 
    ! Ensure a set of three coordinate components is written for 1D, 2D and ND systems
    min_dim = min(size(pos, 1), 3)
    position(1:3) = 0.0_real64
 
    ! xyz-style labels and positions:
    do iatom = 1, n_atoms
      position(1:min_dim) = pos(1:min_dim, iatom)
      write(iunit, '(a10, 3(1x, f15.10))') trim(adjustl(species(iatom))), &
        units_from_atomic(unit_angstrom, position)
    end do
 
    call io_close(iunit)
 
    pop_sub(write_standard_xyz_file)
 
  end subroutine write_standard_xyz_file
 
  subroutine write_canonicalized_xyz_file(dir, fname, space, latt, pos, atoms, box, namespace)
    character(len=*),        intent(in) :: dir
    character(len=*),        intent(in) :: fname
    class(space_t),          intent(in) :: space
    type(lattice_vectors_t), intent(in) :: latt
    real(real64),            intent(in) :: pos(:,:)
    type(atom_t),            intent(in) :: atoms(:)
    class(box_t),            intent(in) :: box
    type(namespace_t),       intent(in) :: namespace
 
    integer                               :: ia
    character(len=:),         allocatable :: info_str, dir_fname
    character(len=LABEL_LEN), allocatable :: species(:)
 
    push_sub(write_canonicalized_xyz_file)
 
    call io_mkdir(dir, namespace)
 
    info_str = trim(space%short_info()) // '; ' // trim(box%short_info(unit_angstrom))
    if (space%is_periodic()) then
      info_str = info_str // '; ' // trim(latt%short_info(unit_angstrom))
    end if
 
    safe_allocate(species(1:size(atoms)))
    do ia = 1, size(atoms)
      species(ia) = atoms(ia)%get_label()
    enddo
 
    dir_fname = trim(dir) // '/' // trim(fname)
    call write_standard_xyz_file(namespace, dir_fname, pos, species, info_str)
 
    safe_deallocate_a(species)
 
    pop_sub(write_canonicalized_xyz_file)
 
  end subroutine write_canonicalized_xyz_file
 
  subroutine write_xsf_geometry_file(dir, fname, space, latt, pos, atoms, mesh, namespace, total_forces)
    character(len=*),        intent(in) :: dir, fname
    class(space_t),          intent(in) :: space
    type(lattice_vectors_t), intent(in) :: latt
    real(real64),            intent(in) :: pos(:,:)
    type(atom_t),            intent(in) :: atoms(:)
    class(mesh_t),           intent(in) :: mesh
    type(namespace_t),       intent(in) :: namespace
    real(real64),  optional,        intent(in) :: total_forces(:,:)
 
    integer :: iunit
    real(real64), allocatable :: forces(:,:)
 
    if (.not. mpi_world%is_root()) return
 
    push_sub(write_xsf_geometry_file)
 
    call io_mkdir(dir, namespace)
    iunit = io_open(trim(dir)//'/'//trim(fname)//'.xsf', namespace, action='write', position='asis')
 
    if (present(total_forces)) then
      safe_allocate(forces(1:space%dim, 1:size(atoms)))
      forces = units_from_atomic(units_out%force, total_forces)
      call write_xsf_geometry(iunit, space, latt, pos, atoms, mesh, forces = forces)
      safe_deallocate_a(forces)
    else
      call write_xsf_geometry(iunit, space, latt, pos, atoms, mesh)
    end if
 
    call io_close(iunit)
 
    pop_sub(write_xsf_geometry_file)
  end subroutine write_xsf_geometry_file
 
  ! ---------------------------------------------------------
  subroutine write_xsf_geometry(iunit, space, latt, pos, atoms, mesh, forces, index)
    integer,                 intent(in) :: iunit
    class(space_t),          intent(in) :: space
    type(lattice_vectors_t), intent(in) :: latt
    real(real64),            intent(in) :: pos(:,:)
    type(atom_t),            intent(in) :: atoms(:)
    class(mesh_t),           intent(in) :: mesh
    real(real64),  optional, intent(in) :: forces(:, :)
    integer,       optional, intent(in) :: index
 
    integer :: idir, idir2, iatom, index_, natoms
    character(len=7) :: index_str
    real(real64) :: offset(space%dim)
    real(real64) :: rlattice(space%dim, space%dim)
 
    push_sub(write_xsf_geometry)
 
    if (present(index)) then
      write(index_str, '(a,i6)') ' ', index
      index_ = index
    else
      write(index_str, '(a)') ''
      index_ = 1
    end if
    natoms = size(pos, dim=2)
 
    ! The corner of the cell is always (0,0,0) to XCrySDen
    ! so the offset is applied to the atomic coordinates.
    ! Along periodic dimensions the offset is -1/2 in reduced coordinates, as
    ! our origin is at the center of the cell instead of being at the edge.
    offset(1:space%dim) = latt%red_to_cart(spread(-m_half, 1, space%dim))
    ! Offset in aperiodic directions:
    do idir = space%periodic_dim + 1, space%dim
      offset(idir) = -(mesh%idx%ll(idir) - 1)/2 * mesh%spacing(idir)
    end do
 
    if (space%is_periodic()) then
      if (index_ == 1) then
        select case (space%periodic_dim)
        case (3)
          write(iunit, '(a)') 'CRYSTAL'
        case (2)
          write(iunit, '(a)') 'SLAB'
        case (1)
          write(iunit, '(a)') 'POLYMER'
        end select
      end if
 
      write(iunit, '(a)') 'PRIMVEC'//trim(index_str)
 
      !Computes the rlattice corresponding to the 3D periodic version of the simulation box
      rlattice = latt%rlattice
      do idir = space%periodic_dim+1, space%dim
        rlattice(:,idir) = rlattice(:,idir)*m_two*mesh%box%bounding_box_l(idir)
      end do
 
      do idir = 1, space%dim
        write(iunit, '(3f12.6)') (units_from_atomic(units_out%length, rlattice(idir2, idir)), idir2 = 1, space%dim)
      end do
 
      write(iunit, '(a)') 'PRIMCOORD'//trim(index_str)
      write(iunit, '(i10, a)') natoms, ' 1'
    else
      write(iunit, '(a)') 'ATOMS'//trim(index_str)
    end if
 
    ! BoxOffset should be considered here
    do iatom = 1, natoms
      write(iunit, '(a10, 3f12.6)', advance='no') trim(atoms(iatom)%label), &
        (units_from_atomic(units_out%length, pos(idir, iatom) - offset(idir)), idir = 1, space%dim)
      if (present(forces)) then
        write(iunit, '(5x, 3f12.6)', advance='no') forces(:, iatom)
      end if
      write(iunit, '()')
    end do
 
    pop_sub(write_xsf_geometry)
  end subroutine write_xsf_geometry
 
! ---------------------------------------------------------
  subroutine write_xsf_geometry_supercell(iunit, space, latt, pos, atoms, mesh, centers, supercell, extra_atom)
    integer,                 intent(in) :: iunit
    class(space_t),          intent(in) :: space
    type(lattice_vectors_t), intent(in) :: latt
    real(real64),            intent(in) :: pos(:,:)
    type(atom_t),            intent(in) :: atoms(:)
    class(mesh_t),           intent(in) :: mesh
    real(real64),            intent(in) :: centers(:, :)
    integer,                 intent(in) :: supercell(:)
    real(real64), optional,         intent(in) :: extra_atom(:)
 
    integer :: idir, idir2, iatom, index_
    character(len=7) :: index_str
    real(real64) :: offset(3)
    integer :: irep, nreplica
 
    push_sub(write_xsf_geometry_supercell)
 
    write(index_str, '(a)') ''
    index_ = 1
 
    nreplica = product(supercell(1:space%dim))
 
    ! The corner of the cell is always (0,0,0) to XCrySDen
    ! so the offset is applied to the atomic coordinates.
    ! Offset in periodic directions:
    offset(1:space%dim) = latt%red_to_cart(spread(-m_half, 1, space%dim))
    offset(1:3) = offset(1:3) + centers(1:3,1)
    ! Offset in aperiodic directions:
    do idir = space%periodic_dim + 1, 3
      offset(idir) = -(mesh%idx%ll(idir) - 1)/2 * mesh%spacing(idir)
    end do
 
    if(space%is_periodic()) then
      if(index_ == 1) then
        select case(space%periodic_dim)
        case(3)
          write(iunit, '(a)') 'CRYSTAL'
        case(2)
          write(iunit, '(a)') 'SLAB'
        case(1)
          write(iunit, '(a)') 'POLYMER'
        end select
      end if
 
      write(iunit, '(a)') 'PRIMVEC'//trim(index_str)
 
      do idir = 1, space%dim
        write(iunit, '(3f12.6)') (units_from_atomic(units_out%length, &
          latt%rlattice(idir2, idir)*supercell(idir)), idir2 = 1, space%dim)
      end do
 
      write(iunit, '(a)') 'PRIMCOORD'//trim(index_str)
      if(.not.present(extra_atom)) then
        write(iunit, '(i10, a)') size(atoms)*nreplica, ' 1'
      else
        write(iunit, '(i10, a)') size(atoms)*nreplica+1, ' 1'
      end if
    else
      write(iunit, '(a)') 'ATOMS'//trim(index_str)
    end if
 
 
    do irep = 1, nreplica
      ! BoxOffset should be considered here
      do iatom = 1, size(atoms)
        write(iunit, '(a10, 3f12.6)', advance='no') trim(atoms(iatom)%label), &
          (units_from_atomic(units_out%length, pos(idir, iatom) + centers(idir, irep) &
          - offset(idir)), idir = 1, space%dim)
        write(iunit, '()')
      end do
    end do
    if(present(extra_atom)) then
      write(iunit, '(a10, 3f12.6)', advance='no') 'X', &
        (units_from_atomic(units_out%length, extra_atom(idir) - offset(idir)), idir = 1, space%dim)
      write(iunit, '()')
    end if
 
    pop_sub(write_xsf_geometry_supercell)
  end subroutine write_xsf_geometry_supercell
 
 
#if defined(HAVE_NETCDF)
  ! ---------------------------------------------------------
  subroutine ncdf_error(func, status, filename, namespace, ierr)
    character(len=*),  intent(in)    :: func
    integer,           intent(in)    :: status
    character(len=*),  intent(in)    :: filename
    type(namespace_t), intent(in)    :: namespace
    integer,           intent(inout) :: ierr
 
    push_sub(ncdf_error)
 
    if (status == nf90_noerr) then
      pop_sub(ncdf_error)
      return
    end if
 
    write(message(1),'(3a)') "NETCDF error in function '" , trim(func) , "'"
    write(message(2),'(3a)') "(reading/writing ", trim(filename) , ")"
    write(message(3), '(6x,a,a)')'Error code = ', trim(nf90_strerror(status))
    call messages_warning(3, namespace=namespace)
    ierr = 5
 
    pop_sub(ncdf_error)
  end subroutine ncdf_error
#endif
 
  ! ---------------------------------------------------------
  subroutine transpose3(in, out)
    real(real64), intent(in)  :: in(:, :, :)
    real(real64), intent(out) :: out(:, :, :)
    integer :: ix, iy, iz
 
    push_sub(transpose3)
 
    do ix = lbound(in, 1), ubound(in, 1)
      do iy = lbound(in, 2), ubound(in, 2)
        do iz = lbound(in, 3), ubound(in, 3)
          out(iz, iy, ix) = in(ix, iy, iz)
        end do
      end do
    end do
 
    pop_sub(transpose3)
  end subroutine transpose3
 
#include "undef.F90"
#include "real.F90"
#include "io_function_inc.F90"
 
#include "undef.F90"
#include "complex.F90"
#include "io_function_inc.F90"
#include "undef.F90"
 
end module io_function_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: