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 on the wiki of octopus. 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 ../../tests/data/methane/exec
/builds/octopus-code/postopus/tests/data/methane/exec
/usr/local/lib/python3.9/dist-packages/IPython/core/magics/osm.py:417: UserWarning: using dhist requires you to install the `pickleshare` library.
  self.shell.db['dhist'] = compress_dhist(dhist)[-100:]
[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/tests/data
/usr/local/lib/python3.9/dist-packages/IPython/core/magics/osm.py:417: UserWarning: using dhist requires you to install the `pickleshare` library.
  self.shell.db['dhist'] = compress_dhist(dhist)[-100:]
[5]:
cd methane
/builds/octopus-code/postopus/tests/data/methane
[6]:
run = Run(".")
[7]:
xa = run.default.scf.density.get_converged("ncdf")

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/tests/data/benzene/exec
/usr/local/lib/python3.9/dist-packages/IPython/core/magics/osm.py:417: UserWarning: using dhist requires you to install the `pickleshare` library.
  self.shell.db['dhist'] = compress_dhist(dhist)[-100:]
[12]:
!head -n 100 parser.log | grep "UnitsOutput"
UnitsOutput = 1
[13]:
cd ..
/builds/octopus-code/postopus/tests/data/benzene
[14]:
runb = Run(".")

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

[15]:
xab = runb.default.scf.density.get_converged(source="ncdf")
xab
[15]:
<xarray.DataArray 'density' (step: 1, x: 47, y: 49, z: 33)>
array([[[[0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         ...,
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.]],

        [[0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         ...,
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.]],

        [[0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         ...,
...
         ...,
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.]],

        [[0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         ...,
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.]],

        [[0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         ...,
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.]]]])
Coordinates:
  * step     (step) int64 1
  * x        (x) float64 -6.9 -6.6 -6.3 -6.0 -5.7 -5.4 ... 5.7 6.0 6.3 6.6 6.9
  * y        (y) float64 -7.2 -6.9 -6.6 -6.3 -6.0 -5.7 ... 6.0 6.3 6.6 6.9 7.2
  * z        (z) float64 -4.8 -4.5 -4.2 -3.9 -3.6 -3.3 ... 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 12.1). 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.get_converged(source="cube")
xa_cube
[19]:
<xarray.DataArray 'density' (step: 1, x: 47, y: 49, z: 33)>
array([[[[0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         ...,
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.]],

        [[0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         ...,
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.]],

        [[0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         ...,
...
         ...,
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.]],

        [[0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         ...,
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.]],

        [[0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         ...,
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.]]]])
Coordinates:
  * step     (step) int64 1
  * x        (x) float64 -13.04 -12.47 -11.91 -11.34 ... 11.34 11.91 12.47 13.04
  * y        (y) float64 -13.61 -13.04 -12.47 -11.91 ... 11.91 12.47 13.04 13.61
  * z        (z) float64 -9.071 -8.504 -7.937 -7.37 ... 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]: