# Output files guidelines {#sec:outputformat}

## General output {#sec:genoutput}

Let\'s say you completed a calculation specified by the following input
file

``` {#code:bn_scf .matlab}
info.calculationType = 'self-consistent'
info.savepath = './results/bn_real_scf'
info.outfile = 'fbn_real_scf.out'
element(1).species = 'B';
element(1).path = './B_AtomicData.mat'
element(2).species = 'N';
element(2).path = './N_AtomicData.mat'
alat = 6.85; % in bohr
atom.xyz = alat*[0 0 0;.25 .25 .25]
atom.element = [1;2]
domain.latvec = alat*[.5 .5 0;.5 0 .5;0 .5 .5]
domain.lowres = 0.2
functional.list = {'XC_GGA_X_PBE','XC_GGA_C_PBE'}
kpoint.gridn = [3 3 3]
kpoint.sigma = 0
mixing.tol = 1e-8*[1 1]
dos.status = 0
option.saveWavefunction = 1
```

The location of the output is determined by the keyword
`info.savepath = ./results/bn_real_scf`. RESCU output can be divided
into two classes:

-   

    **Distributed**

    :   This includes essentially data structures living on some grid
        (describing real space or reciprocal space fields like the
        electron density, the Hartree potential, wavefunctions, etc.)
        and atomic orbital matrices (like the Hamiltonian). This is
        because in large scale systems, any of these quantities is
        potentially too large to fit in the memory of a single node.

-   

    **Global**

    :   This includes pretty much anything else, like the atomic
        coordinates, the lattice vectors, the total energy, etc.

The global quantities are included in the MAT file (in our example
`./results/bn_real_scf.mat`). A MAT-file may be read simply with the
MATLAB command `load`. This is shown below.

``` matlab
>> load('/home/vincentm/Downloads/rescu-2.7.0/examples/results/al_lcao_scf.mat')
>> atom

atom =

struct with fields:

      Ztot: 3
      desc: [1×1 struct]
   element: 1
   nvalence: 3
         xyz: [0 0 0]

>> atom.desc

ans =

struct with fields:

   element: 'Vector containing the element index of every atoms.'
   nvalence: 'atom.nvalence determines the number of valence electrons in self-consistent calculations. The valence is intrinsically linked to the density. By default, atom.nvalence is determined such that the system is charge neutral. It is however possible to choose a value whereby the system is non-neutral.'
         xyz: 'Cartesian coordinates of every atoms in atomic units (Bohr radii).'
```

In each structure like `atom`, there is a `desc` field which contains
descriptions for the parameters and output data.

The distributed quantities are included in the HDF5 file (in our example
`./results/bn_real_scf.h5`). HDF5 files may be parsed with
[h5disp](https://www.mathworks.com/help/matlab/ref/h5disp.html) and read
with [h5read](https://www.mathworks.com/help/matlab/ref/h5read.html) in
MATLAB. For example,

``` 
>> h5disp('results/al_lcao_scf.h5','/rho/rhoOut');
HDF5 al_lcao_scf.h5
Dataset 'rhoOut'
   Size:  729x1
   MaxSize:  729x1
   Datatype:   H5T_IEEE_F64LE (double)
   ChunkSize:  []
   Filters:  none
   FillValue:  0.000000
   Attributes:
      'real':  1.000000
```

## Units

All output quantities are in atomic units, unless otherwise specified.
This means that energies are in Hartree and lengths (area, volumes,
etc.) are in Bohr. For example, one could revise the above
`input file <code:bn_scf>`{.interpreted-text role="ref"} as follows

``` matlab
units.length = 'ang'
alat = 3.6249 % in ang
atom.xyz = alat*[0 0 0;.25 .25 .25]
atom.element = [1;2]
domain.latvec = alat*[.5 .5 0;.5 0 .5;0 .5 .5]
```

The lattice constant is now passed in Angstroms, but the system
specified by the input file is equivalent. In both cases, the output is
in atomic units. For example, one could load and show the lattice
vectors

``` matlab
>> load('/home/vincentm/repos/rescumat/tmp/results/bn_real_scf.mat')
>> domain.latvec

ans =

   3.4250    3.4250         0
   3.4250         0    3.4250
         0    3.4250    3.4250
```

## Output control

As mentioned in section `sec:genoutput`{.interpreted-text role="ref"},
output data can be quite demanding in terms of disk space, and hence not
everything is stored on disk by default. The keywords in table
`tab:saveopts`{.interpreted-text role="numref"} shows various keywords
used to control which output data should be saved

  Keywords                                                                                                              Description                                                                                                                                        HDF5 group
  --------------------------------------------------------------------------------------------------------------------- -------------------------------------------------------------------------------------------------------------------------------------------------- ----------------------------------------------
  option.saveDensity option.saveDFPT option.saveEigensolver option.saveInterpolation option.saveLevel option.saveLCAO   Density (electronic, atomic, etc.) DFPT parameters Eigensolver parameters Interpolation matrices Set to 2 to save everything All LCAO parameters   `rho`
  option.saveLcaoHam option.saveMisc option.saveMixing                                                                  LCAO Hamiltonian Debug option Mixer history                                                                                                        `/LCAO`
  option.savePartialWeights                                                                                             Band decomposition                                                                                                                                 `/LCAO/partialWeight`
  option.savePotential option.saveTau                                                                                   Potential Kinetic energy density                                                                                                                   `/potential`
  option.saveWavefunction                                                                                               Wavefunctions (periodic part)                                                                                                                      `/waveFunction/psi` `/waveFunction/deltaPsi`
  option.saveDensityMatrix                                                                                              Real space density matrix                                                                                                                          `/DensityMatrix/Rcell`
  option.saveDensityMatrixK                                                                                             Density matrix (k-space)                                                                                                                           `/DensityMatrix/kcell`
                                                                                                                        Local DOS                                                                                                                                          `/dos/ldosVal`
  dfpt.saveEPImat                                                                                                       Electron-phonon matrix                                                                                                                             `/dfpt/EPImat`
  dfpt.saveElectronVelocity                                                                                             Electron velocities                                                                                                                                `/dfpt/electronVelocitiy`
  dfpt.saveElectronSE                                                                                                   Electron self-energy (eph scattering)                                                                                                              `/dfpt/electronSE`
  dfpt.saveElectronTau                                                                                                  Electron lifetime (eph scattering)                                                                                                                 `/dfpt/electronTau`
                                                                                                                        Mulliken population matrix                                                                                                                         `/LCAO/mullikenPop`

  : Save keywords

Some quantities like the wavefunctions `/waveFunction/psi` have k-point
and spin indices. The HDF5 datasets are thus numbered consecutively from
1 to N to distinguish them. The first $n_k$ wavefunctions correspond to
the first $n_k$ irreducible k-points specified in `kpoint.ikdirect`. In
collinear spin calculations, the $n_k$ `psi` are for spin up and the
next $n_k$ `psi` datasets are for spin down. The same logic holds for
other quantities like `/LCAO/partialWeight`.

## Grids

In RESCU, there are generally two grids:

-   One for the wavefunctions, the so-called coarse grid with grid
    number given by `domain.cgridn`;
-   One for everything else, the so-called fine grid with grid number
    given by `domain.fgridn`.

In many calculations (e.g. all LCAO calculations), these grids are the
same. But in some cases, they might be different. So when reshaping a
quantity like the density, you should pay attention to use the grid that
matches the number of elements in the array.

Another important point is where the grid is with respect to the unit
cell. In RESCU, the fractional coordinates of grid points along each
dimension are $[du,1]$. Let\'s say you want to work with the fine grid,
you can thus get the 1D grids as follows

``` matlab
nu = domain.fgridn(1);  % # grid points along dimension #1
nv = domain.fgridn(2);  % # grid points along dimension #2
nw = domain.fgridn(3);  % # grid points along dimension #3
u2 = (1:nu)/nu;         % fractional coords of grid along lattice vector #1
v2 = (1:nv)/nv;         % fractional coords of grid along lattice vector #2
w2 = (1:nw)/nw;         % fractional coords of grid along lattice vector #3
```

The full 3D grid can be generated as follows

``` matlab
avec = domain.latvec;
[u2,v2,w2] = ndgrid(u2,v2,w2);
xyz = [u2(:) v2(:) w2(:)]*avec;
x = reshape(xyz(:,1), nu, nv, nw);
y = reshape(xyz(:,2), nu, nv, nw);
z = reshape(xyz(:,3), nu, nv, nw);
```

You can visualize this as follows

``` matlab
figure(); hold on
avec = domain.latvec;
plot_domain_edges(avec);
scatter3(x(:),y(:),z(:), 'filled')
```

and you\'ll see something like figure `grids_al`{.interpreted-text
role="numref"}

<figure>
<img src="data/grids_Al.png" id="grids_al" alt="data/grids_Al.png" />
<figcaption>Grids for an Al calculation. Note the grid starts at <span
class="math inline"><em>d</em><em>x</em></span> not 0.</figcaption>
</figure>

## Arrays

Many quantities are conveniently represented as arrays. For example,
lattice vectors are written as

``` matlab
avec = [[1.,2.,3.];[4.,5.,6.];[7.,8.,9.]]
```

If arrays with more than one dimension are written as vectors, they will
be reshaped automatically. RESCU (like MATLAB) uses the so-called
Fortran order, which means that a vector is reshaped filling the 1st
dimension, then the 2nd, then the 3rd etc. This means that

``` matlab
avec = [1.,2.,3.,4.,5.,6.,7.,8.,9.]
```

is reshaped to

``` matlab
avec = [[1.,4.,7.];[2.,5.,8.];[3.,6.,9.]]
```

which corresponds to a matrix

$$\begin{aligned}
\begin{pmatrix}
1 & 4 & 7\\
2 & 5 & 8\\
3 & 6 & 9\\
\end{pmatrix}
\end{aligned}$$

By contrast, there exist the so-called C ordering where multidimensional
arrays are filled backward, i.e. starting by the last dimension.

This means that

``` matlab
avec = [1.,2.,3.,4.,5.,6.,7.,8.,9.]
```

is reshaped to

``` matlab
avec = [[1.,2.,3.];[4.,5.,6.];[7.,8.,9.]]
```

which corresponds to a matrix

$$\begin{aligned}
\begin{pmatrix}
1 & 2 & 3\\
4 & 5 & 6\\
7 & 8 & 9\\
\end{pmatrix}
\end{aligned}$$

Going back to the example above, suppose you import the density as

``` matlab
>> rho = h5disp('results/al_lcao_scf.h5','/rho/rhoOut');
    Dataset 'rhoOut'
     Size:  15625x1
     MaxSize:  15625x1
     Datatype:   H5T_IEEE_F64LE (double)
     ChunkSize:  []
     Filters:  none
     FillValue:  0.000000
     Attributes:
         'real':  1.000000
```

you should reshape `rho` according to the domain data in
`results/al_lcao_scf.mat`. In the present example,
`domain.fgridn = [9,9,9]` so `rho = reshape(rho, [9,9,9])` will give you
the density as a 3D array.

For the wavefunctions, you would see something like

``` matlab
>> h5disp('results/bn_real_scf.h5','/waveFunction' )
HDF5 bn_real_scf.h5
Group '/waveFunction'
   Dataset 'psi1'
      Size:  15625x12
      MaxSize:  15625x12
      Datatype:   H5T_IEEE_F64LE (double)
      ChunkSize:  []
      Filters:  none
      FillValue:  0.000000
      Attributes:
            'real':  1.000000
   Dataset 'psi2'
      Size:  15625x24
      MaxSize:  15625x24
      Datatype:   H5T_IEEE_F64LE (double)
      ChunkSize:  []
      Filters:  none
      FillValue:  0.000000
      Attributes:
            'real':  0.000000
```

Note that the datasets are not 1D but 2D now, with shape `15625x24`. The
first index runs over the grid points and the second index is for bands.
There are therefore 24 bands of boron nitride in this output.

## Normalization

In RESCU, fields are treated in analogy with analytical fields. This
means that the normalization condition for the density, for instance, is

$$\int_{\Omega} \rho(\mathbf r) d\mathbf r = N_{el}$$

where $\Omega$ is a unit cell and $N_{el}$ is the number of electrons.
Numerically, this translates to (for an aluminum primitive cell)

``` 
>> load('./results/al_lcao_scf.mat')
>> vol = abs(det(domain.latvec));
>> fn = prod(domain.fgridn);
>> dr = vol/fn;
>> rho = h5read('results/al_lcao_scf.h5','/rho/rhoOut');
>> sum(rho)*dr

ans =

   3.0000
```

Similarly, the normalization condition for wavefunctions, or
wavefunction-like quantities like atomic orbitals, is

$$\int_{\Omega} \psi(\mathbf r)^{\dagger} \psi(\mathbf r) d\mathbf r = 1$$

Numerically, this translates to (again for an aluminum primitive cell)

``` 
>> vol = abs(det(domain.latvec));
cn = prod(domain.cgridn);
dr = vol/cn;
psi = h5read('results/al_real_scf.h5','/waveFunction/psi1');
psi'*psi*dr

ans =

Columns 1 through 7

   1.0000   -0.0000   -0.0000   -0.0000    0.0000    0.0000   -0.0000
   -0.0000    1.0000   -0.0000   -0.0000    0.0000    0.0000    0.0000
   -0.0000   -0.0000    1.0000   -0.0000    0.0000    0.0000    0.0000
   -0.0000   -0.0000   -0.0000    1.0000   -0.0000   -0.0000    0.0000
   0.0000    0.0000    0.0000   -0.0000    1.0000    0.0000    0.0000
   0.0000    0.0000    0.0000   -0.0000    0.0000    1.0000   -0.0000
   -0.0000    0.0000    0.0000    0.0000    0.0000   -0.0000    1.0000
   -0.0000   -0.0000   -0.0000   -0.0000    0.0000   -0.0000   -0.0000
   -0.0000    0.0000   -0.0000   -0.0000   -0.0000   -0.0000   -0.0000
   -0.0000    0.0000    0.0000   -0.0000    0.0000   -0.0000    0.0000

Columns 8 through 10

   -0.0000   -0.0000   -0.0000
   -0.0000    0.0000    0.0000
   -0.0000   -0.0000    0.0000
   -0.0000   -0.0000   -0.0000
   0.0000   -0.0000    0.0000
   -0.0000   -0.0000   -0.0000
   -0.0000   -0.0000    0.0000
   1.0000    0.0000   -0.0000
   0.0000    1.0000   -0.0000
   -0.0000   -0.0000    1.0000
```

Kohn-Sham states are orthonormal, which results in the identity matrix.

## Atomic orbitals information

RESCU may solve the Kohn-Sham equation (and others) using linear
combinations of atomic orbitals. Even when using plane waves and real
space grids, it is often useful to project certain quantities on
numerical atomic orbitals to learn about the *character* of, say,
wavefunctions. It is then crucial to know how RESCU stores atomic
orbital information. Take for instance the following system.

``` {caption="Cu3Au.xyz"}
4
Cu3Au
Au 0.0 0.0 0.0
Cu 0.0 0.5 0.5
Cu 0.5 0.0 0.5
Cu 0.5 0.5 0.0
```

``` {caption="Cu3Au_lcao_scf.input"}
LCAO.status = 1
info.calculationType = 'self-consistent'
info.savepath = './results/Cu3Au_lcao_scf'
atom.fracxyz = 'Cu3Au.xyz'
domain.latvec = 3.74*eye(3); units.latvec = 'a'
domain.lowres = 0.250000
element(1).species = 'Au'
element(1).path = './Au_AtomicData.mat'
element(2).species = 'Cu'
element(2).path = './Cu_AtomicData.mat'
functional.list = {'XC_GGA_X_PBE','XC_GGA_C_PBE'}
kpoint.sampling = 'tetrahedron+blochl'
kpoint.gridn = [6 6 6]
LCAO.mulliken = true
```

We set `LCAO.mulliken = true` to save the Mulliken population matrix
upon exit. In the present example, this is a $60 \times 60$ matrix with
one column/row per orbital. How can one know which orbital each column
corresponds to. Loading the output, one sees the structure `LCAO` in the
MATLAB workspace. The orbital data is found in `LCAO.orbInfo`. It is
summarized in `orbinfo`{.interpreted-text role="numref"}.

  Keywords   Description
  ---------- ----------------------------------
  `Aorb`     Atom index
  `Lorb`     Orbital angular momentum
  `Morb`     z-axis orbital angular momentum
  `Oorb`     Orbital index (in PP file)
  `Porb`     Atomic population (neutral atom)
  `Rorb`     Cutoff radius
  `Sorb`     Species index
  `Zorb`     $\zeta$-number

  : Atomic orbital info

Suppose we want to compute the atomic Mulliken population matrix. We
first load the Mulliken population matrix and atom indices as follows.

``` matlab
mul = readH5Array('results/Cu3Au_lcao_scf.h5','/LCAO/mullikenPop1');
aind = LCAO.orbInfo.Aorb
```

Then

``` matlab
amul = zeros(4, 4);
for ii = 1:4
   for jj = 1:4
      amul(ii, jj) = real(sum(sum(mul(aind == ii, aind == jj))));
   end
end
amul
```

``` matlab
amul =

   10.0893    0.3879    0.4371    0.4211
   0.3879   10.0463    0.2450    0.2660
   0.4371    0.2450    9.9917    0.2648
   0.4211    0.2660    0.2648    9.8288
```

Other quantities like the Hamiltonian, partial density of states, etc.
are ordered in the same way.

## Real spherical harmonics

Atomic orbitals are defined as the product of a radial basis function
and a spherical harmonic as follows

$$\phi_{I\mu}(\mathbf{r}) = \zeta_{\mu}(|\mathbf{r}-\mathbf{R}_I|)Y_{l_\mu m_\mu}(\widehat{\mathbf{r}-\mathbf{R}_I})$$

There are several conventions for spherical harmonics. In what follows
we\'ll use

$$Y_\ell^m( \theta , \varphi ) = (-1)^m \sqrt{\frac{(2\ell+1)}{4\pi}\frac{(\ell-m)!}{(\ell+m)!}} \, P_{\ell}^m ( \cos{\theta} ) \, e^{i m \varphi }$$

For performance reasons, RESCU uses real spherical harmonics instead of
the usual complex spherical harmonics. Real spherical harmonics are
defined in `tab:realsh`{.interpreted-text role="numref"} below.

  $l$   $m$                                                             $\tilde Y_l^m$   
  ----- ----- --------------------------------------------------------- ---------------- --------------------------------------------------------------------------------------
  0     0     $Y_0^0$                                                   $=$              $\frac{1}{2} \sqrt{\frac{1}{\pi}}$
  1     -1    $i \sqrt{\frac{1}{2}} \left( Y_1^{- 1} + Y_1^1 \right)$   $=$              $\sqrt{\frac{3}{4 \pi}} \cdot \frac{y}{r}$
  1     0     $Y_1^0$                                                   $=$              $\sqrt{\frac{3}{4 \pi}} \cdot \frac{z}{r}$
  1     1     $\sqrt{\frac{1}{2}} \left( Y_1^{- 1} - Y_1^1 \right)$     $=$              $\sqrt{\frac{3}{4 \pi}} \cdot \frac{x}{r}$
  2     -2    $i \sqrt{\frac{1}{2}} \left( Y_2^{- 2} - Y_2^2\right)$    $=$              $\frac{1}{2} \sqrt{\frac{15}{\pi}} \cdot \frac{x y}{r^2}$
  2     -1    $i \sqrt{\frac{1}{2}} \left( Y_2^{- 1} + Y_2^1 \right)$   $=$              $\frac{1}{2} \sqrt{\frac{15}{\pi}} \cdot \frac{y z}{r^2}$
  2     0     $Y_2^0$                                                   $=$              $\frac{1}{4} \sqrt{\frac{5}{\pi}} \cdot \frac{- x^2 - y^2 + 2 z^2}{r^2}$
  2     1     $\sqrt{\frac{1}{2}} \left( Y_2^{- 1} - Y_2^1 \right)$     $=$              $\frac{1}{2} \sqrt{\frac{15}{\pi}} \cdot \frac{z x}{r^2}$
  2     2     $\sqrt{\frac{1}{2}} \left( Y_2^{- 2} + Y_2^2 \right)$     $=$              $\frac{1}{4} \sqrt{\frac{15}{\pi}} \cdot \frac{x^2 - y^2 }{r^2}$
  3     -3    $i \sqrt{\frac{1}{2}} \left( Y_3^{- 3} + Y_3^3 \right)$   $=$              $\frac{1}{4} \sqrt{\frac{35}{2 \pi}} \cdot \frac{\left( 3 x^2 - y^2 \right) y}{r^3}$
  3     -2    $i \sqrt{\frac{1}{2}} \left( Y_3^{- 2} - Y_3^2 \right)$   $=$              $\frac{1}{2} \sqrt{\frac{105}{\pi}} \cdot \frac{xy z}{r^3}$
  3     -1    $i \sqrt{\frac{1}{2}} \left( Y_3^{- 1} + Y_3^1 \right)$   $=$              $\frac{1}{4} \sqrt{\frac{21}{2 \pi}} \cdot \frac{y (4 z^2 - x^2 - y^2)}{r^3}$
  3     0     $Y_3^0$                                                   $=$              $\frac{1}{4} \sqrt{\frac{7}{\pi}} \cdot \frac{z (2 z^2 - 3 x^2 - 3 y^2)}{r^3}$
  3     1     $\sqrt{\frac{1}{2}} \left( Y_3^{- 1} - Y_3^1 \right)$     $=$              $\frac{1}{4} \sqrt{\frac{21}{2 \pi}} \cdot \frac{x (4 z^2 - x^2 - y^2)}{r^3}$
  3     2     $\sqrt{\frac{1}{2}} \left( Y_3^{- 2} + Y_3^2 \right)$     $=$              $\frac{1}{4} \sqrt{\frac{105}{\pi}} \cdot \frac{\left( x^2 - y^2 \right) z}{r^3}$
  3     3     $\sqrt{\frac{1}{2}} \left( Y_3^{- 3} - Y_3^3 \right)$     $=$              $\frac{1}{4} \sqrt{\frac{35}{2 \pi}} \cdot \frac{\left( x^2 - 3 y^2 \right) x}{r^3}$

  : Real spherical harmonics definitions
