# How can I efficiently write formatted data to a Gaussian Cube file?

I am working with data files that contain information about electronic wave functions. The data format is specified by the quantum chemistry program Gaussian. I need to import the data, manipulate it, then output it so it can be read by other programs. The file format is specified here, and I put an example file here. The problem is that outputting the data to a file takes an unconscionably long time, and I want to know if that time can be reduced or if it is as fast as it can be.

In the example below, the data file is small (less than 2mb), but I am currently working with much larger files, where just writing the file takes over a minute and ends up being the rate-limiting step in my calculations.

The data can be read in via

    header = Import["http://pastebin.com/raw.php?i=JXtLghmD", "Text"];
cubedata = ToExpression[Import["https://gist.github.com/anonymous/f37baaa89e8a17d39f12/raw/d872df53f9251827ef4a7ec6f00087dd7f94552d/form_MO_28"]];


If I only need to use the data inside of Mathematica, then the easiest and fastest thing to do is to write it to a binary format.

    Export["cube_28.nc", {"Cube" -> cubedata,"HeaderText" -> header}] // AbsoluteTiming


When I run the above code, it takes two thirds of a second to execute and the resulting file is about 1mb. But if I need to write the data so molecular visualization software can read it, I have to write a text file. Since there is no good printf in Mathematica, this was improvised:

    OutForm[num_?NumericQ, width_Integer, ndig_Integer, OptionsPattern[]] :=
Module[{mant, exp, val}, {mant, exp} = MantissaExponent[num];
mant = ToString[NumberForm[mant, {ndig, ndig}]];
exp = If[Sign[exp] == -1, "-", "+"] <> IntegerString[exp, 10, 2];
val = mant <> "E" <> exp;

Module[{moltxt, nx, ny, nz},
moltxt = OpenWrite[cubeFileName, FormatType -> FortranForm];
WriteString[moltxt, "\n"];
{nx, ny, nz} = Dimensions[cubeData];
Do[Do[Do[WriteString[moltxt, OutForm[cubeData[[x, y, z]], 13, 5]];
If[Mod[z, 6] == 0, WriteString[moltxt, "\n"]];, {z, 1, nz}];
WriteString[moltxt, "\n"];, {y, 1, ny}];, {x, 1, nx}];
Close[moltxt];];


And now I can write the data to a formatted Gaussian cube file via

    WriteCube["cube_28.cube", header, cubedata] // AbsoluteTiming


which takes eight seconds to execute and results in a 1.68mb file.

Does creating the formatted ASCII output have to take twelve times as long to create? Is there any way OutForm and WriteCube could be written more efficiently?

• Does your molecular visualization program accept binary files and if yes can you find the format? Apr 16, 2013 at 1:07
• Would that it did, I'd be home free. But I don't believe it does. I use VMD, or Avogadro to visualize the data and they both accept cube files as input. Apr 16, 2013 at 1:10
• I see, for VMD you have to alter the plugin to read binary files too. I dont know how feasible that might be. It's amazing the only binary format supported by OpenBabel is cdx! I think optimizing your WriteCube command is the best bet, at least to reduce the time. ASCII files will always be many times bigger than binary you can't do anything about that. Apr 16, 2013 at 1:30
• I don't have time to test it right now but I expect better performance by forming larger blocks of data to pass to WriteString, rather than using it inside three nested Do loops. Apr 16, 2013 at 1:34

This provides a modest speed-up, and if I didn't mess up it should be equivalent:

WriteCube2[cubeFileName_?StringQ, headerTxt_?StringQ, cubeData_] :=
Module[{stream},
stream = OpenWrite[cubeFileName, FormatType -> FortranForm];
Map[
WriteString[stream, ##, "\n"] & @@
Riffle[OutForm[#, 13, 5] & /@ #, "\n", {7, -1, 7}] &,
cubeData, {2}
];
Close[stream];
]


Theoretically it would be more efficient to pass even larger blocks of data to WriteString (that is, call it fewer times), but using this method reveals that with this code the vast majority of the time is spent in the OutForm operation so additional attention should be focused there. Vector/array optimizations may be possible by writing OutForm to apply to the entire tensor at once.

• Yes, I will take a look at OutForm to see what I can do there. There is an implementation of printf lurking somewhere in the stack overflow archives, and it was much much slower than the OutForm function. I tried replacing OutForm[#,13,5] with ScientificForm[#, NumberFormat -> (Row[{#1, "e", #3, " "}] &)] but the result isn't correct. I want it to print the same number of characters regardless of how many significant digits it thinks the number has. Apr 16, 2013 at 15:44
• Hi, I wanted to re-use your code, but it looks like a magic. Could you please explain where all these numerical constants, like 13, 5, 7, -1, 7 pop up from? Feb 27, 2021 at 12:26
• @yarchik I haven't looked at the code in years. 13 and 5 are parameters of Jason's code. {7, -1, 7} is a Riffle specification to insert every seventh element, including a trailing element. Compare Riffle[Range[12], X, {7, -1, 7}] to Riffle[Range[12], X, 7]. The data format specifies six columns of data, and inserting "\n" in this way appears to have been my method of line-breaking. Mar 6, 2021 at 18:52