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

Question

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;
       StringJoin@PadLeft[Characters[val], width, " "]];

    WriteCube[cubeFileName_?StringQ, headerTxt_?StringQ, cubeData_] := 
      Module[{moltxt, nx, ny, nz}, 
       moltxt = OpenWrite[cubeFileName, FormatType -> FortranForm];
       WriteString[moltxt, headerTxt];
       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?

Answer 1

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];
   WriteString[stream, headerTxt, "\n"];
   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.