# Generating numpy-compatible electron density from Molden / QC calculation

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

# 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?

• 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 Jun 8 at 2:54
• 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 Jun 8 at 12:21
• @Wesley - maybe you have some time now? I'd love to see an answer, so I'm giving a bounty. :-) Jun 23 at 23:07
• @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
– Tyberius
Jun 24 at 17:42
• 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.
– Tyberius
Jul 6 at 19:03

# ORBKIT

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:

ok.init()
# 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.

• I'll also mention that there's a command-line program, e.g. orbkit -i file.molden.gz --otype=cb --adjust_grid=5 0.2 -p 4 which will write a Gaussian-style cube file, which can be read into numpy or pandas easily Jun 25 at 13:33