Detailed overview#

Introduction#

Postopus is a post-processing tool for Octopus (POSTprocessing for OctoPUS). It provides a user-friendly interface to find and read data written by Octopus throughout a simulation and offers common operations to evaluate this data.

Octopus takes an input file describing the systems and the simulation parameters. The input file is called inp and is a text file with no extension. The command octopus would then have to be executed in the folder that contains this inp file. Since octopus doesn’t allow custom names for the input file, a possible project structure could look like the following:

├── benzene                    # Project folder
│   ├── benzene.xyz            # Geometry file or other supporting files
│   └── inp                    # Input file
├── h-atom
│   └── inp
├── he
│   └── inp
├── methane                     # Project folder in case of a multi-stage calculation
│   ├── calculation_gs          # Ground state calculation
│   │   └── inp
│   ├── calculation_td          # Time dependent calculation
│   │   └── inp
│   └── inp                     # Input file for the whole calculation (The other files must be placed here one by one in each stage)
└── recipe
    └── inp

To then run one of these simulations one would run the command octopus in the root of the respective folder.

Running the simulation#

As mentioned before, running a simulation involves two steps:

Change the directory to the project folder (that contains the inp file)
Run the command octopus(optionally store the octopus output in a log file by calling octopus > out_gs.log 2>&1).

The above two steps could theoretically be executed in a separate shell. However, since the Jupyter notebook has terminal capabilities, it is recommended to execute all steps of a workflow (from defining the input file to the final analysis) in a single notebook (when computationally feasible), as we will do in the following. The notebook then serves as a record of all the steps taken to achieve a particular result. This increases the reproducibility of your work and makes it easier for others to understand and reuse your conclusions. For the latter reason, we also recommend using !pip freeze in the notebook to keep track of the software versions of the used Python packages and to print out the Octopus version used to generate the data with octopus -v.

[1]:

!octopus -v

octopus 16.0 (git commit 7e864b450a)

[2]:

!pip freeze

accessible-pygments==0.0.5
alabaster==1.0.0
anyio==4.9.0
argon2-cffi==23.1.0
argon2-cffi-bindings==21.2.0
arrow==1.3.0
ase==3.25.0
asttokens==3.0.0
async-lru==2.0.5
attrs==25.3.0
babel==2.17.0
beautifulsoup4==4.13.4
bleach==6.2.0
bokeh==3.7.2
certifi==2025.1.31
cffi==1.17.1
cftime==1.6.4.post1
charset-normalizer==3.4.1
click==8.1.8
cloudpickle==3.1.1
colorcet==3.1.0
comm==0.2.2
contourpy==1.3.2
cycler==0.12.1
dask==2025.3.0
debugpy==1.8.14
decorator==5.2.1
defusedxml==0.7.1
docutils==0.21.2
executing==2.2.0
fastjsonschema==2.21.1
fonttools==4.57.0
fqdn==1.5.1
fsspec==2025.3.2
h11==0.14.0
holoviews==1.20.2
httpcore==1.0.8
httpx==0.28.1
idna==3.10
imagesize==1.4.1
importlib_metadata==8.6.1
ipykernel==6.29.5
ipython==9.1.0
ipython_pygments_lexers==1.1.1
ipywidgets==8.1.6
isoduration==20.11.0
jedi==0.19.2
Jinja2==3.1.6
json5==0.12.0
jsonpointer==3.0.0
jsonschema==4.23.0
jsonschema-specifications==2024.10.1
jupyter==1.1.1
jupyter-console==6.6.3
jupyter-events==0.12.0
jupyter-lsp==2.2.5
jupyter_client==8.6.3
jupyter_core==5.7.2
jupyter_server==2.15.0
jupyter_server_terminals==0.5.3
jupyterlab==4.4.0
jupyterlab_pygments==0.3.0
jupyterlab_server==2.27.3
jupyterlab_widgets==3.0.14
kiwisolver==1.4.8
linkify-it-py==2.0.3
locket==1.0.0
Markdown==3.8
markdown-it-py==3.0.0
MarkupSafe==3.0.2
matplotlib==3.10.1
matplotlib-inline==0.1.7
mdit-py-plugins==0.4.2
mdurl==0.1.2
mistune==3.1.3
narwhals==1.35.0
nbclient==0.10.2
nbconvert==7.16.6
nbformat==5.10.4
nbsphinx==0.9.7
nest-asyncio==1.6.0
netCDF4==1.7.2
notebook==7.4.0
notebook_shim==0.2.4
numpy==2.2.5
overrides==7.7.0
packaging==25.0
pandas==2.2.3
pandocfilters==1.5.1
panel==1.6.2
param==2.2.0
parso==0.8.4
partd==1.4.2
pexpect==4.9.0
pickleshare==0.7.5
pillow==11.2.1
platformdirs==4.3.7
pooch==1.8.2
-e git+https://gitlab-ci-token:glcbt-66_z2qz92VWzov64xZCAqUv@gitlab.com/octopus-code/postopus.git@3c5bdf1fd47489f4a65e61254d5a448966144c91#egg=postopus
prettytable==3.16.0
prometheus_client==0.21.1
prompt_toolkit==3.0.51
psutil==7.0.0
ptyprocess==0.7.0
pure_eval==0.2.3
pyarrow==19.0.1
pycparser==2.22
pydata-sphinx-theme==0.16.1
Pygments==2.19.1
pyparsing==3.2.3
python-dateutil==2.9.0.post0
python-json-logger==3.3.0
pytz==2025.2
pyvista==0.45.0
pyviz_comms==3.0.4
PyYAML==6.0.2
pyzmq==26.4.0
referencing==0.36.2
requests==2.32.3
rfc3339-validator==0.1.4
rfc3986-validator==0.1.1
rpds-py==0.24.0
scipy==1.15.2
scooby==0.10.0
Send2Trash==1.8.3
six==1.17.0
sniffio==1.3.1
snowballstemmer==2.2.0
soupsieve==2.7
Sphinx==8.1.3
sphinxcontrib-applehelp==2.0.0
sphinxcontrib-devhelp==2.0.0
sphinxcontrib-htmlhelp==2.1.0
sphinxcontrib-jsmath==1.0.1
sphinxcontrib-qthelp==2.0.0
sphinxcontrib-serializinghtml==2.0.0
stack-data==0.6.3
terminado==0.18.1
tinycss2==1.4.0
toolz==1.0.0
tornado==6.4.2
tqdm==4.67.1
traitlets==5.14.3
types-python-dateutil==2.9.0.20241206
typing_extensions==4.13.2
tzdata==2025.2
uc-micro-py==1.0.3
uri-template==1.3.0
urllib3==2.4.0
vtk==9.4.2
wcwidth==0.2.13
webcolors==24.11.1
webencodings==0.5.1
websocket-client==1.8.0
widgetsnbextension==4.0.14
xarray==2025.3.1
xrft==1.0.1
xyzservices==2025.1.0
zipp==3.21.0

[3]:

!mkdir -p "examples/interference"

[4]:

cd examples/interference/

/builds/octopus-code/postopus/docs/notebooks/examples/interference

We now create the input file (inp). To do this we use the magic command %%writefile inp of the Jupyter Notebook to write the contents of the cell to a file called inp (our input file) in the current directory.

[5]:

%%writefile inp

stdout = "td_stdout.txt"
stderr = "td_stderr.txt"

CalculationMode = td
ExperimentalFeatures = yes
FromScratch = yes

%Systems
  'Maxwell' | maxwell
%

Maxwell.ParStates = no

# Maxwell box variables
lsize_mx = 10.0

Maxwell.BoxShape = parallelepiped

%Maxwell.Lsize
 lsize_mx | lsize_mx | lsize_mx
%

dx_mx = 0.5

%Maxwell.Spacing
 dx_mx | dx_mx | dx_mx
%

# Maxwell calculation variables
%MaxwellBoundaryConditions
 plane_waves | plane_waves | plane_waves
%

%MaxwellAbsorbingBoundaries
 not_absorbing | not_absorbing | not_absorbing
%

# Output variables
OutputFormat = plane_x + plane_y + plane_z + axis_x + axis_y + axis_z

# Maxwell output variables
%MaxwellOutput
 electric_field
 magnetic_field
 maxwell_energy_density
 trans_electric_field
%
MaxwellOutputInterval = 10
MaxwellTDOutput = maxwell_energy + maxwell_total_e_field + maxwell_total_b_field

# Time step variables
TDSystemPropagator = prop_expmid
dt = 1 / ( sqrt(c^2/dx_mx^2 + c^2/dx_mx^2 + c^2/dx_mx^2) )
TDTimeStep = dt
TDPropagationTime = 0.35

# laser propagates in x direction
k_1_x  =  0.707107
k_1_y  = -0.707107
k_2_x  = -0.447214
k_2_y  = -0.223607
E_1_z  =  0.5
E_2_z  =  0.5
pw_1   =  5.0
pw_2   =  7.5
ps_1_x = -sqrt(1/2) * 20.0
ps_1_y =  sqrt(1/2) * 20.0
ps_2_x =  sqrt(2/3) * 20.0
ps_2_y =  sqrt(1/3) * 20.0

%MaxwellIncidentWaves
  plane_wave_mx_function | electric_field | 0 | 0 | E_1_z | "plane_waves_function_1"
  plane_wave_mx_function | electric_field | 0 | 0 | E_2_z | "plane_waves_function_2"
%

%MaxwellFunctions
  "plane_waves_function_1" | mxf_cosinoidal_wave | k_1_x | k_1_y | 0 | ps_1_x | ps_1_y | 0 | pw_1
  "plane_waves_function_2" | mxf_cosinoidal_wave | k_2_x | k_2_y | 0 | ps_2_x | ps_2_y | 0 | pw_2
%

Writing inp

Assuming you have octopus in your PATH:

[6]:

!octopus

Loading Data with Postopus#

To load data with Postopus the path to the output directory of the Octopus simulation is required. In this folder, all output data, as well as the input file inp are expected. Data is found automatically and can be discovered by the user by listing all found systems/fields/etc or using auto-completion at run time, e. g. when using Jupyter Notebook.

The entry point for users to Postopus is the Run class.

[7]:

from postopus import Run

run = Run()

The Run object discovers available data on the file system and builds a data structure allowing access. If no path to a directory containing the inp file is passed to the run object (i.e. run = Run()), then the run object searches in the current working directory. The run object can also be instantiated with a specific path (i.e. run = Run("path/to/inpfile_directory/")). In general, the data structure allows choosing data with the following syntax:

run.systemname.calculationmode.output_name

Parameters set in italics must be replaced with values that mostly correspond to values set in the input file. A closer look at those will be taken in the following sections.

System selection#

The first parameter to select is the system’s name. Octopus allows to simulate multiple systems at once, as well as so-called “multisystem”s which build a hierarchy of encapsulated systems.

Checking out the “Systems” block in the inp, the Maxwell system can be found:

%Systems
  'Maxwell' | maxwell
%

One system with the name “Maxwell” of type “maxwell”. The types here are relevant for Octopus, for us the system names are of interest.

Be aware that simulation with Octopus is also possible without setting any systems. In that case, the system’s type will be set (by Octopus) to “electronic_system”. As Postopus requires a name for this system, it will be named “default” (while not having a name in Octopus). Also, the “default” system will always exist, as it is used to store global parameters read from the inp, but will never contain any data when the “Systems” block is defined in inp.

Besides reading these names from the inp, it also is possible to access this via Postopus. Use:

[8]:

run

[8]:

Run('/builds/octopus-code/postopus/docs/notebooks/examples/interference'):
Found systems:
    'Maxwell'
    'default'
Found calculation modes:
    'td'

System data - Calculation modes and subsytems#

To load data from a system, we now call: run.Maxwell. This gives a list of the calculation modes which are available via Postopus:

[9]:

run.Maxwell

[9]:

System(name='Maxwell', rootpath='/builds/octopus-code/postopus/docs/notebooks/examples/interference'):
Found calculation modes:
    'td'

This is expected, as the CalculationMode variable in the inp is set to “td”. As the time-dependent calculation (“td”) required a previous self-consistent field simulation (“scf”) this data also could be present in the output folder. If this would be the case, one could see

System(name='Maxwell', rootpath='.'):
Found calculation modes:
    'scf'
    'td'

as output and select between these two. For multisystem examples like the celestial_bodies tutorial, we would also have the keys Moon, Earth and Sun as subsystems, for example run.SolarSystem.Earth and run.SolarSystem.Moon or, if three levels of nesting are used, run.SolarSystem.Earth.Terra and run.SolarSystem.Earth.Luna.

Outputs#

Getting a list of all available outputs can be done with:

[10]:

run.Maxwell.td

[10]:

CalculationModes('td'):
Found 7 outputs:
    b_field
    e_field
    e_field_trans
    maxwell_energy
    maxwell_energy_density
    total_b_field
    total_e_field

[11]:

run.Maxwell.td.maxwell_energy

[11]:

StaticOutput(name='maxwell_energy')

Call the output to get the provided data

[12]:

run.Maxwell.td.maxwell_energy()

[12]:

	t	Total	E**2	B**2	Total+Boundaries	Boundaries	Transversal	Longitudinal	Incident Waves
Iter
0	0.000000	0.000000e+00	0.000000e+00	0.000000e+00	0.000000	0.000000	0.0	0.0	0.000000
1	0.002107	0.000000e+00	0.000000e+00	0.000000e+00	0.033354	0.033354	0.0	0.0	0.016677
2	0.004213	8.108214e-09	6.026031e-09	2.082183e-09	0.061151	0.061151	0.0	0.0	0.030575
3	0.006320	7.915800e-08	3.900977e-08	4.014822e-08	0.101315	0.101315	0.0	0.0	0.050657
4	0.008426	6.407517e-07	1.436090e-07	4.971427e-07	0.155032	0.155031	0.0	0.0	0.077516
...	...	...	...	...	...	...	...	...	...
162	0.341263	2.216257e-03	1.188263e-03	1.027994e-03	0.002216	0.000000	0.0	0.0	0.001028
163	0.343370	2.162107e-03	1.165330e-03	9.967770e-04	0.002162	0.000000	0.0	0.0	0.000997
164	0.345477	2.106947e-03	1.125132e-03	9.818150e-04	0.002107	0.000000	0.0	0.0	0.000982
165	0.347583	2.049632e-03	1.072135e-03	9.774969e-04	0.002050	0.000000	0.0	0.0	0.000977
166	0.349690	1.989371e-03	1.012270e-03	9.771001e-04	0.001989	0.000000	0.0	0.0	0.000977

167 rows × 9 columns

Octopus produces output files accross multiple folders on the file system. Postopus tries to group them together and provide them in a single object.

As example, Octopus outputs the files “e_field-x”, “e_field-y” and “e_field-z” (with multiple extensions “.x=0”, “.y=0”, “.x=0,y=0”, …) for each n step of the simulation in “output_iter/td.0000000”, “output_iter/td.0000010”, “output_iter/td.0000020”, …

In Postopus all these files are united and provided by the “e_field” output:

[13]:

run.Maxwell.td.e_field

[13]:

FieldOutput(name='e_field'):
Available sources:
    '.x=0'
    '.x=0,y=0'
    '.x=0,z=0'
    '.y=0'
    '.y=0,z=0'
    '.z=0'

[14]:

data = run.Maxwell.td.e_field(source=".x=0")
data

The data is provided as xarray.Dataset here. By selecting data.sel(t=0.21, method="nearest") one could access the files in “”output_iter/td.0000010”. By accessing the component data.vy one would read from e_field-y.x=0.

Note that the data is accessed “lazily”. This means, as long one does not work with the data no files are accessed. Only when the data is needed (e.g. when doing data.values, data.min(), …) loading the files will be invoked. Also only the files which match to the selection are used. In data.sel(t=0.21, method="nearest").values files in “output_iter/td.0000000” won’t be touched. Be aware that accessing the values before selecting specific steps may take a while when dealing with long simulations.

The data provided by postopus is in most cases either a pandas DataFrame (as in run.Maxwell.td.maxwell_energy()) or an xarray DataArray/Dataset (as in run.Maxwell.td.e_field(source=".x=0")). Outputs returning an xarray object are referred to as “field”.

For outputs where different sources (file extensions) are available the source has to be provided as parameter.

Plotting#

Plotting the data is well integrated with common packages, such as matplotlib and holoviews:

[26]:

import matplotlib.pyplot as plt

e_field_plane = run.Maxwell.td.e_field(source="z=0").sel(t=0.15, method="nearest")

plt.imshow(e_field_plane.vz.values);

../_images/notebooks_Detailed_Overview_47_0.png

Xarray also provides a convenient .plot method. It uses matplotlib internally and makes use of the additional metadata such as dimensions and coordinates to create axis labels and a colorbar. For more details please refer to the plotting Documentation.

[27]:

e_field_plane.vz.plot(x="x");

../_images/notebooks_Detailed_Overview_49_0.png

To generate dynamic and/or higher dimensional plots we recommend using holoviews, details can be found in the holoviews tutorial:

[28]:

import holoviews as hv

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

[29]:

# Note for web users: You should have an active notebook session to interact with the plot
e_field_over_time = run.Maxwell.td.e_field(source="z=0")
hv_ds = hv.Dataset(e_field_over_time.vz)
hv_im = hv_ds.to(hv.Image, kdims=["x", "y"])
hv_im

[29]:

Postprocessing#

To manipulate Xarray data is in general fairly simple, as it integrates many built-in methods from the scipy and NumPy libraries among others (see the documentation), having the advantage that the operations that involve spatial manipulation are more intuitive (see also for example the xrft tutorial). The results of the computations are itself xarrays, so it is still possible to do use all the plotting methods presented above.

[30]:

integrated_field = e_field_over_time.vz.integrate(coord="y")

[31]:

integrated_field

As we integrated one of the coordinates, we are dealing now with 1D data that evolves in time, thus we are not going to use holoviews.Images, but holoviews.Curves.

[32]:

# Note for web users: You should have an active notebook session to interact with the plot
hv_ds_int = hv.Dataset(integrated_field)
hv_im_int = hv_ds_int.to(hv.Curve, kdims=["x"], dynamic=True)
hv_im_int

[32]: