jjgoings / two-component_cubes

Python script to prepare complex, two-component Gaussian cube files to be read by GaussView or Pymol, etc.

Cubes

Basic use

This is a python class to read and write out Gaussian cube files, with a particular emphasisis on extracting complex and two-component molecular orbitals.

Example:

Say you have a complex 2c MO cube called twoc.cube. If you have cube.py in the same directory, you could write a script like

from cube import Cube
twoc = Cube('twoc.cube')
twoc.write_out('twoc_ra.cube',data='RA')
twoc.write_out('twoc_ia.cube',data='IA')
twoc.write_out('twoc_rb.cube',data='RB')
twoc.write_out('twoc_ib.cube',data='IB')

Or, if you just cloned the whole directory you can run cube.py as-is

$ python cube.py

(You can modify the behavior of cube.py underneath the if__name__ == '__main__()':)

This would create an object called twoc, which you could then dump out new cube files that are readable by GaussView corresponding to the real alpha (data='RA'), real beta (data='RB'), imaginary alpha (data='IA'), and imaginary beta (data='IB'). The default for write_out is real alpha.

The code should be smart enough to distinguish between MO and density cube files.

Usually there is no problem using GaussView for RHF and UHF anyway, so the code isn't yet generalized to all cases, but I'll work that out when I have time.

Vector quantities: Getting magnetization densities and electrostatic field gradients

If you generated your GHF cubes using the VSPIN (in the development version only of cubegen, I believe), you can also write out magnetization densities.

So, for example you ran the following:

cubegen 0 vspin ghf-vspin-example.fchk ghf-vspin-example.cube 40 h

Then you could do, e.g.

magnetization = Cube('ghf-vspin-example.cube',vspin=True)
magnetization.write_out('Mz.cube',data='Mz')

This would dump the z-component of the magnetization density to a new cube called Mz.cube. You can do this for the norm of the magnetization data='N', as well as the x and y components of the magnetization (data='Mx' and data='My', respectively). Densitites are always real.

There is also a new function to plot some vector quantities, if you have Mayavi installed. These are handled through the plot_property method. So continuing the example above,

magnetization = Cube('ghf-vspin-example.cube',vspin=True)
magnetization.plot_property('MAG')

would plot the magnetization vector at different points in 3D space. It's in beta, so you may need to fiddle around with the source to get it to do what you want.

Similarly, if you have an electrostatic potential generated by, for example,

cubegen 0 potential=scf esp-example.fchk esp-example.cube 80 h

Then you can do

ESP = Cube('esp-example.cube')
ESP.plot_property('ESP')

which will plot the gradient (a vector) of your electrostatic potential. This is helpful for visualizing the electric field generated by a molecule with a dipole moment, for example.