I'd like to use a variety of Python tools (e.g., matplotlib and Mayavi) to analyze orbitals plot electron density contours, etc. Most of these expect to work on numpy matrices.

My ideal workflow would be something like this:

import packageX

molecule = packageX.read('file.molden')
# specify grid options
density = packageX.calculateDensity(molecule)

Are there Python packages that make it easy to calculate electron density, spin density, and orbitals on a NumPy grid?

  • 4
    $\begingroup$ An old colleague of mine had a nice way to do it with a package he called Lightspeed. There was an additional program needed to read the molden orbitals, and then it called light speed to do some density calculation. I assume you can do the same with PSI4 $\endgroup$
    – Cody Aldaz
    Commented Jun 8, 2021 at 2:54
  • 7
    $\begingroup$ iodata + gbasis modules of HORTON can do this. I can't make the script for you right now, but, HORTON has your back on this one. github.com/theochem/horton $\endgroup$
    – Wesley
    Commented Jun 8, 2021 at 12:21
  • 2
    $\begingroup$ @Wesley If its any additional incentive, this question is also one of our top 10 unanswered questions, for which I have a standing 50 point bounty posted on Meta $\endgroup$
    – Tyberius
    Commented Jun 24, 2021 at 17:42
  • $\begingroup$ Since you wound up finding an answer yourself, you would be eligible for the bounty I mentioned to Wesley. It may lead to an alternative answer to the question or it may just make up some of the loss from your bounty. Just wanted to check with you before putting the bounty on this post. $\endgroup$
    – Tyberius
    Commented Jul 6, 2021 at 19:03

1 Answer 1



This uses the ORBKIT package.

Gunter Hermann, Vincent Pohl, Jean Christophe Tremblay, Beate Paulus, Hans-Christian Hege, and Axel Schild, "ORBKIT: A Modular Python Toolbox for Cross-Platform Postprocessing of Quantum Chemical Wavefunction Data", J. Comput. Chem. 2016, 37, 1511-1520.

import orbkit as ok

ok.options.filename = '/path/to/86.molden' # also directly reads file.molden.gz if needed
ok.options.itype = 'molden'
ok.options.no_output = True # we'll directly get a Numpy array without generating a file
ok.options.adjust_grid = [5, 0.1] # 5 bohr padding on all sides, 0.1 bohr resolution
ok.options.numproc = 4 # or whatever you want - uses Python subprocess

data = ok.run_orbkit() # prints useful log bits in Jupyter

data.shape # returns (186, 179, 282) for my example file

If you want to do this in a loop, you'll need to set the ok.options every time:

# Set some options (again)
ok.options.filename = '/path/to/new/file.molden'
ok.options.itype = 'molden'
ok.options.no_output = True
ok.options.adjust_grid = [5, 0.1] # 5 bohr padding on all sides, 0.1 bohr resolution
ok.options.numproc = 4 # or whatever you want - uses Python subprocess
# run orbkit again
data_2 = ok.run_orbkit()

Orbkit easily handles the following formats:

  • Molden
  • wfn/wfx
  • fchk
  • Interface to cclib parser (for other formats)

Beyond electron density it's pretty easy to generate:

  • orbitals
  • spin density
  • derivatives (analytical first and second order)
  • Laplacian

More can be found on the ORBKIT website

I had to patch ORBKIT to read the Molden files I had (https://github.com/orbkit/orbkit/pull/39) but you can install via pip install orbkit or python setup.py install once you download.


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .