Visualize Orbitals

import exatomic
from exatomic.base import resource         # Easy access to static files
from exatomic import UniverseWidget as UW  # The visualization system

Quantum chemistry codes don’t always make it easy to get all the necessary information but provided the outputs have the data and the parser is implemented to handle the print-out from a specific code, this is what the API looks like for orbital viewing. - First you parse the output - must contain full basis set specification and a matrix of basis function coefficients

ed = exatomic.myqmcode.Output('/path/to/my/output')

• Inspect the parsed data. Does it look correct?

ed.basis_set         # houses primitive exponents, coefficients
ed.momatrix          # contains the basis function coefficients
ed.basis_set_order   # indices specifying full basis set order

• ed is an Editor, an object that makes parsing files fun!

  • ed must be converted to a universe for the magic to happen

uni = ed.to_universe()

• A Universe has an add_molecular_orbitals method

uni.add_molecular_orbitals?   # to see all keyword arguments
uni.add_molecular_orbitals()  # tries to guess smart defaults

• Don’t forget python uses 0-based indexing, so check your vector values!

Gaussian API

from exatomic import gaussian

uni = gaussian.Output(resource('g09-ch3nh2-631g.out')).to_universe()
fni = gaussian.Fchk(resource('g09-ch3nh2-631g.fchk')).to_universe()

uni.atom # Inspect the geometry

	set	Z	x	y	z	frame	symbol
atom
0	0	6	0.100257	1.346217	0.000000	0	C
1	1	7	0.100257	-1.431666	0.000000	0	N
2	2	1	-1.847313	2.016818	0.000000	0	H
3	2	1	1.074043	2.016820	1.686645	0	H
4	2	1	1.074043	2.016820	-1.686645	0	H
5	2	1	-0.802058	-2.053050	-1.562856	0	H
6	2	1	-0.802058	-2.053050	1.562856	0	H

uni.basis_set_order.head() # And the first few basis functions

	center	L	ml	shell	frame
chi
0	0	0	0	0	0
1	0	0	0	1	0
2	0	1	1	3	0
3	0	1	-1	3	0
4	0	1	0	3	0

# Let's add the first 10 molecular orbitals
uni.add_molecular_orbitals(vector=range(10))
fni.add_molecular_orbitals(vector=range(10))

Evaluating 28 basis functions once.
Timing: compute orbitals -     2.81s.
Evaluating 28 basis functions once.
Timing: compute orbitals -     0.37s.

The UniverseWidget accepts up to 9 universes

• Play around with the buttons:

  • Clear removes contents from the scene. (Press Fill to get the geometry back)

  • Active Scenes controls which scenes are controlled by the buttons (default all)

  • Image allows for saving PNGs of the scenes (can specify a directory and file names)

  • Camera allows to link the camera between scenes (load and save are buggy)

  • Fill switches between WebGL fragment shaders and fancy three.js Spheres

  • Axis adds a unit axis (often hidden within an atom if centered at the origin)

  • Animate will play dynamics (if a universe has multiple frames)

  • Fields expands to show:

    • Isosurfaces tab shows the molecular orbitals

    • Contours tab shows contour lines in x, y, or z axes

• Be sure to use the close button to clean up the javascript side of things

UW(uni, fni)

Recursive solid harmonics means no hard-coded maximum angular momentum

uo2 = gaussian.Output(resource('g09-uo2.out')).to_universe()
uo2.add_molecular_orbitals()
UW(uo2)

Evaluating 141 basis functions once.
Timing: compute orbitals -     1.06s.

Molcas API

• You must include the BSSHOW keyword in the SEWARD module or GATEWAY with a higher than normal print threshold to print the basis set specification

• Molcas prints out Orb files which contain basis function coefficients so we use 2 Editors

from exatomic import molcas

mol = molcas.Output(resource('mol-uo2-anomb.out'))
orb = molcas.Orb(resource('mol-uo2-anomb.scforb'))
# Just attach it to the universe
# mol.momatrix = orb.momatrix
# mol.orbital = orb.orbital

# Alternatively there's a convenience method on molcas.Output
mol.add_orb(resource('mol-uo2-anomb.scforb'))  # adds momatrix and orbital
uni = mol.to_universe()

uni.add_molecular_orbitals(vector=range(40, 60))

Evaluating 69 basis functions once.
Timing: compute orbitals -     3.37s.

UW(uni)

NWChem API

from exatomic import nwchem

nw = nwchem.Output(resource('nw-ch3nh2-631g.out')).to_universe()

nw.add_molecular_orbitals(vector=range(20))
UW(nw)

Evaluating 28 basis functions once.
Timing: compute orbitals -     0.43s.

• If you are interested in helping, that’s great!

• The Editor objects deal with file IO and contain parse_methods specific to each dataframe.

• In order to add_molecular_orbitals, the Editor must parse the Atom, BasisSet, BasisSetOrder and MOMatrix dataframes.

• The requirements for these dataframes are specified in their docstrings.

• The actual interface with the widget occurs on the Universe object, so the Editor must be converted to a Universe.

• Therefore, it makes sense to subclass the exatomic.Editor to gain access to the .to_universe() method.

• Keep parse_dataframe methods modular if possible to isolate problems when they arise.