export CUBE files

[jrussell]

Dear ATK users, Please pardon my ignorance, as I am a new user! I would like to export electron density in a CUBE format. I can do this in VNL just fine, and I can also have the data outputted into my output file. Here's what I do in my job file to ask for the CUBE of electron density:

Code

electron_density = ElectronDensity(
    configuration=bulk_configuration,
    spin=Spin.Sum,
    )
nlsave('density.nc', electron_density)
nlprint(electron_density)

I am having a problem though -- the coordinate information for the atoms is not included in the cubefile and I have to add it manually. Is there a script, or a way to have to get the cube file in the proper format? Here's the problem:

Code

CUBE FILE.
OUTER LOOP: X, MIDDLE LOOP: Y, INNER LOOP: Z
    0 0.000000e+00 0.000000e+00 0.000000e+00
   28 0.000000e+00 1.558687e-01 1.558687e-01
   28 1.558687e-01 0.000000e+00 1.558687e-01
   28 1.558687e-01 1.558687e-01 0.000000e+00

  1.04643e-04  6.94395e-01  1.03695e+00  8.03578e-01  5.12214e-01  3.04839e-01
......(the rest of the file)

In the CUBE FILE, It is saying that there are 0 atoms in my unit cell on line 3. I have to fix this manually, and then manually insert the missing atomic number & coordinates on line 7. Any ideas on how to get this without the extra editing? Thanks for your help, John

[zh]

Please see here:
http://quantumwise.com/forum/index.php?topic=1310.0
http://quantumwise.com/forum/index.php?topic=447.0
http://quantumwise.com/forum/index.php?topic=1216.0

[jrussell]

Please see here: http://quantumwise.com/forum/index.php?topic=1310.0 http://quantumwise.com/forum/index.php?topic=447.0 http://quantumwise.com/forum/index.php?topic=1216.0

Hi Zh, thanks for your reply. I am still stuck! I try the solution for this one: http://quantumwise.com/forum/index.php?topic=1310.0 This is using cube_rev_eng.py script. I get an error:

Code

[jrussell@login1 test]$ atkpython cube_rev_eng.py vnl.cube > out.cube
Traceback (most recent call last):
  File "cube_rev_eng.py", line 20, in <module>
    header_data[0].append(int(line[0:5].strip()))
ValueError: invalid literal for int() with base 10: ''

I don't have any idea why this didn't work. I tried the next solution for this one: http://quantumwise.com/forum/index.php?topic=447.0 This is using export_cube.py included with the second forum answer. I had to modify the script that calls export_cube.py because ATK.TwoProbe is no longer valid. Also restoreSelfConsistentCalculation() is no longer valid. And also calculateElectronDensity() is no longer valid.

Code

#!/bin/env atkpython
from NanoLanguage import *
from export_cube import *

#scf = restoreSelfConsistentCalculation('checkpoint.nc')
configuration = nlread("checkpoint.nc")[0]
electron_density = ElectronDensity(
    configuration=configuration,
    spin=Spin.Sum,
    )
export_cube_file(electron_density,'out.cube')

And I got this result:

Code

Traceback (most recent call last):
  File "dump_cube.py", line 11, in <module>
    export_cube_file(electron_density,'out.cube')
  File "/home/jrussell/atk/li2o/LiO/test/export_cube.py", line 37, in export_cube_file
    origin = grid._cell().origin()
AttributeError: 'ElectronDensity' object has no attribute '_cell'

Perhaps this is an old solution no longer valid for ATK 12.8.2? This old stuff is very misleading. The last forum link didn't seem to apply. Thanks very much for your help, and any advice is greatly appreciated! Best, John

[Anders Blom]

Those links above would not be relevant for your problem. I'll provide a solution for you in a bit, just need to study the cube file format a bit. What are you importing into?

[jrussell]

Hi Anders, thanks for your help! I am importing the file into the Bader analysis program by G. Henkelman, et. al.
(found here: http://theory.cm.utexas.edu/bader/ ). It needs a Gaussian CUBE file format as an input
(file format here: http://paulbourke.net/dataformats/cube/ ). I am trying to calculate the atomic charges from the
electron density (mulliken charges are too far off for my application).

I also use VMD and Chemcraft visualisation software to render isosurfaces generated from CUBE files for publications,
so it is very useful to be able to export into the CUBE file format.

Thanks again for your help,
John

[jrussell]

Hi Anders,

I also have an extra question -- is it possible to increase the density of grid points in the CUBE file?
For example, it would be helpful to be able to get grid spacing of  12 points/Bohr. I didn't see how
to change the grid density in the ElectronDensity keyword.

Thanks again,
John

[Anders Blom]

Here is a first attempt to solve both your issues. I made a simple script with all functions, but if you plan to use this a lot you may want to split it up into a module you import.

Anyway, first step is to see if the produced cube file is understood by the program you import into. I hope the script is clear, note that you only need to understand and modify the last few lines ("User part"), the rest is generic.

To compute what mesh cut off you need for a particular point density, use the formula at http://quantumwise.com/documents/manuals/latest/ReferenceManual/index.html/ref.numericalaccuracyparameters.html.

Have fun!

[jrussell]

Hi Anders, Thanks for the awesome script! This is very helpful! I had to make some minor tweaks --

Code

    header = """CUBE FILE. DATA VALUES IN UNITS OF %s, LENGTHS IN BOHR
OUTER LOOP: X, MIDDLE LOOP: Y, INNER LOOP: Z
%5i%18.6e%18.6e%18.6e
%5i%18.6e%18.6e%18.6e
%5i%18.6e%18.6e%18.6e
%5i%18.6e%18.6e%18.6e
"""

Here I changed the number format from %13.6 to %18.6e to prevent numbers from running together. The second tweak was to add a column to the atom section. There are a total five columns: atomic number, charge, x, y, and z. Here I just added the fifth column and set the charge to zero as I am using the electron density to get Bader charges anyways.

Code

    if number_of_atoms:
        coordinates = gridvalues._configuration().cartesianCoordinates().inUnitsOf(Bohr)
        elements = gridvalues._configuration().elements()
        zero = 0.0
        for e,c in zip(elements,coordinates):
            stream.write("%5i%4.1f%18.6e%18.6e%18.6e\n" % (e.atomicNumber(),zero,c[0],c[1],c[2]))
    else:
        stream.write("\n")

But that was really minor stuff that I could figure out This is really useful, thanks again for your help! Best, John

[Anders Blom]

Good work! I missed the charge column, thanks for fixing that.

We implemented the column widths according to the original CUBE file spec, and if you imagine reading the file in a Fortran code hard-coded to the spec, the fact that the numbers float together is less of any issue but the number of characters per field is. But you'll quickly find out and it's easy to change the export code as you did.