Gaussian includes a standalone utility for generating cubes from the data in a formatted checkpoint file (equivalent to the previous Cube keyword). The utility is named cubegen, and it has the following syntax:
cubegen nprocs kind fchkfile cubefile npts format cubefile2
The parameters, which are not case-sensitive, have the following meanings:
Number of shared memory processors used for electrostatic potential calculations. A value of 0 is equivalent to 1 (it is the default). Note that this parameter must be included if other parameters are specified. Previously this parameter was used to specify the amount of memory to allocate. The GAUSS_MEMDEF environment variable should be used instead.
A keyword specifying the type of cube to generate:
|MO=n||Molecular orbital n. The keywords Homo, Lumo, All, OccA (all alpha occupied), OccB (all beta occupied), Valence (all valence orbitals) and Virtuals (all virtual orbitals) may also be used in place of a specific orbital number. There is no default for n, and an error will occur if it is omitted. AMO and BMO can be similarly used to select only alpha or beta orbitals (respectively). For open shell systems, Homo selects both alpha and beta orbitals.|
|Density=type||Total density of the specified type. The type keyword is one of the single density selection options that are valid with the Density keyword: SCF, MP2, CC, CI, and so on (note that Current is not supported). The fdensity, falpha and fbeta forms request the use of full instead of frozen core densities. The default is SCF.|
|Spin=type||Spin density (difference between α and β densities) of the specified type.|
|Alpha=type||Alpha spin density of the specified type.|
|Beta=type||Beta spin density of the specified type.|
|Potential=type||Electrostatic potential using the density of the specified type.|
|Gradient||Compute the density and gradient.|
|Laplacian||Compute the Laplacian of the density (∇2ρ).|
|NormGradient||Compute the norm of the density gradient at each point.|
|CurrentDensity=I||Magnitude of the magnetically-induced (GIAO) current density, where I is the applied magnetic field direction (X, Y or Z).|
|ShieldingDensity=IJN||Magnetic shielding density. I is the direction of the applied magnetic field, J is the direction of the induced field (X, Y or Z), and N is the number of the nucleus for which the shielding density (GIAO) is to be calculated.|
Name of the formatted checkpoint file. cubegen will prompt for this filename if it is not specified.
Name of the output cube file; test.cube is the default if it is not explicitly specified (i.e., specifying the name of the checkpoint file does not change the default cube filename).
Number of points per side in the cube. A value of 0 selects the default value of 803 points distributed evenly over a rectangular grid generated automatically by the program (not necessarily a cube). Positive values of npts similarly specify the number of points per “side”; e.g., 100 specified a grid of 1,000,000 (1003) points.
The values -2, -3 and -4 correspond to the keywords Coarse, Medium and Fine and to values of 3 points/Bohr, 6 points/Bohr and 12 points/Bohr (respectively). Negative values of npts ≤-5 specify spacing of npts*10-3 Angstroms between points in the grid.
A value of -1 says to read the cube specification from the input stream or from a second cube file (see below), according to the following format:
IFlag, X0, Y0, Z0 Output unit number and initial point. N1, X1, Y1, Z1 Number of points and step-size in the 1-direction. N2, X2, Y2, Z2 Number of points and step-size in the 2-direction. N3, X3, Y3, Z3 Number of points and step-size in the 3-direction.
IFlag is the output unit number. If IFlag is less than 0, then a formatted file will be produced; otherwise, an unformatted file will be written.
If N1<0 the input cube coordinates are assumed to be in Bohr, otherwise, they are interpreted as Angstroms. |N1| is used as the number of X-direction points in any case; N2 and N3 specify the number of points in the Y and Z directions, respectively. Note that the three axes are used exactly as specified; they are not orthogonalized, so the grid need not be rectangular.
The value -5 says to read in an arbitrary list of points from standard input. If you enter this input by hand, terminate the input with an end-of-file (i.e., Ctrl-D under Unix). Alternatively, you can redirect standard input to a file containing the list of points (do not place a blank line or Ctrl-D at the end of the file). The coordinates of the points must be in Angstroms.
Format of formatted output files: h means include header (this is the default); n means don't include header. This parameter is ignored when unformatted cube files are produced.
If specified, the size for the generated cube is taken from this file. This option is useful when creating cubes for later arithmetic operations such as difference densities.
In order for the cube dimensions to be taken from the specified cube file, the npts parameters must be -1, and the specified file must have been created with a header.
If no parameters are specified, cubegen will prompt for fchkfile and run the following:
cubegen 1 density=scf test.cube 80 h
This generates a file named test.cube (with header) containing the SCF density in a rectangular grid of 803 points.
Output File Formats
In general, all values in the cube file are in atomic units, regardless of the input units. The exception to this are the cube files created with a user-provided arbitrary list of points (npts=-5). In this case, the output cube will contain one line per point, with the first 3 values of each line being the Cartesian coordinates of the point in Angstroms (just as given in the input), followed by the values of the reported properties (density, ESP, etc) in atomic units.
Densities, the norm of the density gradient, the Laplacian of the density, and the potential are scalars (i.e. one value per point, NVal=1). A gradient cube contains the density plus the vector gradient of the density, so it has four values per point (NVal=4): i.e. the value of the density plus the X, Y, and Z components of its gradient.
Cube files have one row per record (i.e., N1*N2 records each of length N3*NVal). For formatted output, each row is written out in format (6E13.5). In this case, if N3*NVal is not a multiple of six, then there may be blank space in some lines.
For example, schematically, cube files will look like this:
|NAtoms, X-Origin, Y-Origin, Z-Origin, NVal||NVal is the # of values/point|
|N1, X1, Y1, Z1||# of increments in the slowest running direction|
|N2, X2, Y2, Z2|
|N3, X3, Y3, Z3||# of increments in the fastest running direction|
|IA1, Chg1, X1, Y1, Z1||Atomic number, charge, and coordinates of the first atom|
|IAn, Chgn, Xn, Yn, Zn||Atomic number, charge, and coordinates of the last atom|
|(N1*N2) records, each length N3*NVal||Values of the density at each point in the grid|
Note that a separate write is used for each record.
For molecular orbital output, NAtoms will be less than zero, and an additional record follows the data for the final atom (in format 10I5 if the file is formatted):
|NMO, (MO(I),I=1,NMO)||Number of MOs and their numbers|
If NMO orbitals were evaluated, then each record is NMO*N3 long and has the values for all orbitals at each point together (in this case NMO=NVal).
Reading Cube Files with Fortran Programs
If one wishes to read the values of a cube file back into an array dimensioned X(NVal,N3,N2,N1), code like the following Fortran loop may be used:
Do 10 I1 = 1, N1 Do 10 I2 = 1, N2 Read(n,'(6E13.5)') ((X(IVal,I3,I2,I1),IVal=1,NVal)I3=1,N3) 10 Continue
where n is the unit number corresponding to the cube file.
Last updated on: 03 Aug 2023. [G16 Rev. C.01]