Units#

Internally, octopus uses atomic units for computing values. Nonetheless, the user can choose the output units with the inp parameter UnitsOutput. The default is atomic [units], but one can use ev_angstrom as well. In the latter case, the spatial dimensions will be measured in Angstrom and the energies will be measured in eV. In the atomic case, the spatial dimensions will be measured in bohr, while all the other variables will be measured in atomic units. More on this in the octopus documentation. The election of the inp parameter of the user is reflected in the parser.log after the run is finished. We retrieve the unit information from there.

We will compare the benzene and methane test examples.

Note: We will assume that you already know from other tutorials that defining the input file and running octopus in the same notebook is recommended. We will not do it here, we just assume the data exsists in the right folder.

Methane example (atomic)#

For the methane example, we’ve left the default atomic units, internally Octopus interprets it as a 0:

[1]:

cd ../octopus_data/methane/exec

/builds/octopus-code/postopus/docs/octopus_data/methane/exec

[2]:

!head -n 100 parser.log | grep "UnitsOutput"

UnitsOutput = 0         # default

We load the data:

[3]:

from postopus import Run
import holoviews as hv
from holoviews import opts  # For setting defaults

hv.extension("bokeh", "matplotlib")  # Allow for interactive plots

[4]:

cd ../../

/builds/octopus-code/postopus/docs/octopus_data

[5]:

cd methane

/builds/octopus-code/postopus/docs/octopus_data/methane

[6]:

run = Run(".")

[7]:

xa = run.default.scf.density("ncdf").isel(step=-1)

Here, we can see that the xarray has the units au. In the following we will show how to access these units and also the units of the individual coordinates:

[8]:

xa.units

[8]:

'au'

[9]:

xa.y.units

[9]:

'Bohr'

If we use holoviews for plotting, the coordinates are automatically labeled according to the units of the xarray. Unfortunately, labeling the color bar needs to be done manually. For more information about plotting with holoviews, look at the [holoviews]((holoviews_with_postopus.ipynb) tutorial.

[10]:

hv_ds = hv.Dataset(xa)
hv_im = hv_ds.to(hv.Image, kdims=["x", "y"], dynamic=True)
hv_im.opts(
    colorbar=True,
    width=500,
    height=400,
    cmap="seismic",
    clabel=f"{xa.name} ({xa.units})",
)

[10]:

Benzene example (ev_angstrom)#

For the benzene example, we specified UnitsOutput = ev_angstrom in the inp file. This is interpreted as the number 1 internally in Octopus:

[11]:

cd ../benzene/exec

/builds/octopus-code/postopus/docs/octopus_data/benzene/exec

[12]:

!head -n 100 parser.log | grep "UnitsOutput"

UnitsOutput = 1

[13]:

cd ..

/builds/octopus-code/postopus/docs/octopus_data/benzene

[14]:

runb = Run(".")

We will do the same exploration as we did before with methane:

[15]:

xab = runb.default.scf.density(source="ncdf").isel(step=-1)
xab

[15]:

<xarray.DataArray 'density' (x: 47, y: 49, z: 33)> Size: 608kB
[75999 values with dtype=float64]
Coordinates:
    step     int64 8B 16
  * x        (x) float32 188B -6.9 -6.6 -6.3 -6.0 -5.7 ... 5.7 6.0 6.3 6.6 6.9
  * y        (y) float32 196B -7.2 -6.9 -6.6 -6.3 -6.0 ... 6.0 6.3 6.6 6.9 7.2
  * z        (z) float32 132B -4.8 -4.5 -4.2 -3.9 -3.6 ... 3.6 3.9 4.2 4.5 4.8
Attributes:
    units:    eV_Ångstrom
[16]:

xab.units

[16]:

'eV_Ångstrom'

[17]:

xab.x.units

[17]:

'Ångstrom'

[18]:

hv_dsb = hv.Dataset(xab)
hv_imb = hv_dsb.to(hv.Image, kdims=["x", "y"])
hv_imb.opts(
    colorbar=True,
    width=500,
    height=400,
    cmap="seismic",
    clabel=f"{xab.name} ({xab.units})",
)

[18]:

Cube format as the exception#

The only (known) exception to this rule is the cube format. The format itself specifies that the values should be calculated in atomic units (if the number of voxels in a dimension is positive, which to the best of our knowledge, is always the case in Octopus 13.0). Thus, the cube output will always be in bohr, independently of what the UnitsOutput parameter says. Actually, xcrysden has also an analogous specification but is ignored by octopus for now. So, we will read xcrysden with the unit that the user specifies in the inp, but this may change in the future

[19]:

xa_cube = runb.default.scf.density(source="cube")
xa_cube

[19]:

<xarray.DataArray 'density' (step: 4, x: 47, y: 49, z: 33)> Size: 2MB
[303996 values with dtype=float64]
Coordinates:
  * step     (step) int64 32B 5 10 15 16
  * x        (x) float64 376B -13.04 -12.47 -11.91 -11.34 ... 11.91 12.47 13.04
  * y        (y) float64 392B -13.61 -13.04 -12.47 -11.91 ... 12.47 13.04 13.61
  * z        (z) float64 264B -9.071 -8.504 -7.937 -7.37 ... 7.937 8.504 9.071
Attributes:
    units:    au
[20]:

xa_cube.units

[20]:

'au'

[21]:

xa_cube.x.units

[21]:

'Bohr'

We will clearly see that the scaling is different in comparison to the last plot, due to the change of units:

[22]:

hv_ds_cube = hv.Dataset(xa_cube)
hv_im_cube = hv_ds_cube.to(hv.Image, kdims=["x", "y"])
hv_im_cube.opts(
    colorbar=True,
    width=500,
    height=400,
    cmap="seismic",
    clabel=f"{xa_cube.name} ({xa_cube.units})",
)

[22]: