/*
                       Copyright (c) 1995 by:
        Leif Laaksonen , Centre for Scientific Computing , ESPOO, FINLAND
            Confidential unpublished property of Leif Laaksonen
                        All rights reserved


      This is a program to convert the output using the "cube"
      command from the Gaussian program to a plot format recognized
      by gOpenMol.




      Run this program on the output from GaussianXX using the
      'cube' keyword in GaussianXX.

      This program produces a coordinate file and a plot file, 
      that can be read into SCARECROW or gOpenMol.

      The running of this program is quite obvious (se later in this
      file).

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 This program converts the Gaussian output using the 'cube'
 command to a form understandable to gopenmol.
 Usage:
 gcube2plt -iinput.file -ooutput.file
 Options:  -mXXX , where XXX is the molecular orbital number to be
                   placed in the plot file
           -p      prevent the output of the coordinate file
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


      If you need any help please feel free to contact:

      Leif.Laaksonen@csc.fi


+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Taken from the;
                 Gaussian 94 User's Reference Manual
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

> Cube keyword


DESCRIPTION

The Cube properties keyword can be used to evaluate molecular orbitals, the 
electrostatic potential, the electron density, density gradient, the norm of 
the density gradient, and Laplacian of the density over a 3 dimensional grid 
(cube) of points. By default, Cube evaluates the electron density (corresponding to
the Density option). Which density is used is controlled by the Density keyword; 
use Density=Current to evaluate the cube over the density from a correlated or 
Cl-Singles wavefunction instead of the default Hartree-Fock density.

Note that only one of the available quantities can be evaluated within any one job 
step. Save the checkpoint file (using %Chk), and include Guess=(Reod,Only) 
Density=Checkpoint in the route section of a subsequent job (or job step) in 
order to evaluate a different quantity without repeating any of the other steps
 of the calculation.

Gaussian 94 provides reasonable defaults for grids, so Cube no longer requires 
that the cube be specified by the user. However, the output filename must still
always be provided (see below).

Alternatively, Cube may be given a parameter specifying the number of points to
use per "side" (the default is 80). For example, Cube=100 specifies a grid of 
1,000,000 points ( 100 ), evenly distributed over the rectangular grid generated
by the program (which is not necessarily a cube). In addition, the input
format used by earlier versions of Gaussian is still supported; Cube=Cards
indicates that a grid will be input. It may be used to specify a grid of
arbitrary size and shape.

The files create by Cube can be manipulated using the cubman utility, described
in chapter 5.

Note that Pop=None will inhibit cube file creation.

INPUT FORMAT

When the user elects to provide it, the grid information is read from the input
stream. The first line-required for all Cube jobs-gives a file name for the cube
file. Subsequent lines, which are included only with Cube=Cards, must conform to
format (15,3F12.6), according to the following syntax:

Output-file-name              Required in all Cube jobs.
IFlag, X0, Y0, Z0             Output unit number and initial point.
N1, X1, Y1, Z1                Number of points and step-size in the X-direction.
N2, X2, Y2, Z2                Number of points and step-size in the Y-direction.
N3, X3, Y3, Z3                Number of points and step-size in the Z-direction.

If IFlag is positive, the output file is unformatted; if it is negative, 
the output file is formatted. If N10, the output is unformatted. If
IFlag<0, the output is formatted. All values in the cube file are in atomic units,
regardless of the input units.

For density and potential grids, unformatted files have one row per record 
(i.e., N1 * N2 records each of length N3). For formatted output, each row is 
written out in format (6E13.5). In this case, if N3 is not a multiple of six,
then there may be blank space in some lines.

The norm of the density gradient is also a scalar (i.e., one value per point),
and is written out in the same manner. Density+gradient grids are similar, but
with two writes for each row (of lengths N3 and 3*N3). Density+gradient+Laplacian
grids have 3 writes per row (of lengths N3, 3*N3, and N3).

For example, for a density cube, the output file looks like this:

NAtoms, X-Origin, Y-Origin, Z-Origin
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, Chgl, X1, Y1, Z1       Atomic number, charge, and coordinates of thefirst atom
...
IAn, Chgn, Xn, Yn, Zn       Atomic number, charge, and coordinates of the last atom

(N1 * N2) records, each of length N3  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 lOI5 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.

READING CUBE FILES WITH FORTRAN PROGRAMS

If one wishes to read the values of the density or potential back into an 
array dimensioned X(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(I3,I2,I1),I3=1,N3)
10    Continue

where n is the unit number corresponding to the cube file.

If the origin is (X0,Y0,Z0), and the increments (X1,Y1,Z1), then point 
(I1,I2,I3) has the coordinates:

X-co0rdinate: X0+(I1-1)*X1+(I2-1)*X2+(I3-1)* X3
Y-coordinate: Y0+(I1-1)*Y1+(I2-1)*Y2+(I3-1)* Y3
Z-coordinate: Z0+(I1-1)*Z1+(I2-1)*z2+(I3-1)* Z3

The output is similar if the gradient or gradient and Laplacian of the charge 
density are also requested, except that in these cases there are two or three 
records, respectively, written for each pair of I1, I2 values. Thus, if the 
density, gradient, and Laplacian are to be read into arrays D(N3,N2,N1),
G(3,N3,N2,N1), RL(N3,N2,N1) from a formatted output file, a correct set of 
Fortran loops would be:

      Do 10 I1 = 1, N1
      Do 10 I2 = 1, N2
        Read(n,'(6F13.5)') (D(I3,I2,I1),I3=1,N3)
        Read(n,'(6F13.5)') ((G(IXYZ,I3,I2,I1),IXYZ=1,3), I3=1,N3)
        Read(n,'(6F13.5)') (RL(I3,I2,I1),I3=1,N3)
10    Continue

where again n is the unit number corresponding to the cube file.


OPTIONS

Density      Compute just the density values. This is the default.

Potential    Compute the electrostatic potential at each point.

Gradient     Compute the density and gradient.

Laplacian    Compute the density, gradient, and Laplacian ofthe density 
             ((nabla)2(rho)).

NormGradient Compute the norm of the density gradient at each point.

Orbitals     Compute the values of one or more molecular orbitals at each point.

FrozenCore   Remove the SCF core density. This is the default for the density,
             and is not allowed for the potential.

Full         Evaluate the density including all electrons.

Total        Use the total density. This is the default.

Alpha        Use only the alpha spin density.

Beta         Use only the beta spin density.

Spin         Use the spin density (difference between alpha and beta densities).

Cards        Read grid specification from the input stream (as described above).

Example 1: Produce the total electron density.

x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-
#p rhf/6-31g* 5d test geom=modela  cube=(density,read) FormCheck=OptCart

Gaussian Test Job 257:
Density cube

0,1
o h f

test257.cube
  -51        -2.0        -2.0        -1.0
   40         0.1         0.0         0.0
   40         0.0         0.1         0.0
   20         0.0         0.0         0.1
	   
x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-


Example 2: Produce the data to plot molecular orbitals

x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-
#p rhf/6-31g* 5d test geom=modela  cube=(orbitals) FormCheck=OptCart

Gaussian Test Job 257:
Density cube

0,1
o h f

leif257.cube
  -51        -2.0        -2.0        -1.0
   40         0.1         0.0         0.0
   40         0.0         0.1         0.0
   20         0.0         0.0         0.1
ALL

x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-


!!!!!!!!!!!!!!!!!!O B S E R V E!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
The cube data has to be given as an orthogonal x-, y-, z-coordinate system
in such a way that the x-axis comes first and the z-axis is given as the 
last one. This means that x is the slowest running coordinate and z is the
fastes running coordinate.

This program produces also a coordinate file in the CHARMM 'crd' format
which can be read by gOpenMol or SCARECROW.

!!!!!!!!!!!!!!!!!!O B S E R V E!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


*/

#include 
#include 
#include 
#include 
#include 

#define GAUSSIAN_TYPE   200
#define MAX_TITLE_LINES   5
#define Rabs(a)    ( ( a ) > 0.0 ? (a) : -(a))
#define SMALL 1.e-05
#define BUFF_LEN 256
#define BOHR_RADIUS 0.52917715  /* conversion constant */  

#define FWRITE(value_p , size)    { Items = \
                                 fwrite(value_p, size , 1L , Output_p);\
                                 if(Items < 1) {\
                     printf("?ERROR - in writing contour file (*)\n");\
                     return(1);}}

#define FWRITEN(value_p , num , size) { Items = \
                                fwrite(value_p, size , num , Output_p);\
                   if(Items < 1) {\
                     printf("?ERROR - in writing contour file (**)\n");\
                     return(1);}}

/* ................................................................... */
char *Usage = 
"This program converts the Gaussian output using the 'cube'\n\
 command to a form understandable to gopenmol.\n\
 Usage:\n\
 gcube2plt -iinput.file -ooutput.file\n\
 Options:  -mXXX , where XXX is the molecular orbital number to be\n\
                   placed in the plot file\n\
           -p      prevent the output of the coordinate file\n\
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n\
 The 'Cube' input for GaussianXX has to be defined as an orthogonal\n\
 coordinate system, where the x coordinate is the slowest running and\n\
 and z the fastest running coordinate\n";
/* ................................................................... */

   char *AtomSymbols={"\
Ac  Ag  Al  Am  As  Au  B   Ba  Be  Bi  Br  C   Ca  Cd  \
Ce  Cl  Co  Cr  Cs  Cu  D   Dy  Er  Eu  F   Fe  Ga  Gd  \
Ge  H   Hf  Hg  Ho  I   In  Ir  K   La  Li  Lu  Mg  Mn  \
Mo  N   Na  Nb  Nd  Ni  Np  O   Os  P   Pa  Pb  Pd  Pm  \
Po  Pr  Pt  Pu  Ra  Rb  Re  Rh  Ru  S   Sb  Sc  Se  Si  \
Sm  Sn  Sr  Ta  Tb  Tc  Te  Th  Ti  Tl  Tm  U   V   W   \
Y   Yb  Zn  Zr  He  Ne  Ar  Kr  Xe  Rn  "};


    int AtomSymbol_p[] =
       { 89 , 47 , 13 , 95 , 33 , 79 , 5  , 56 , 4  , 83 , 35 , 6  , 20 , 48 ,
         58 , 17 , 27 , 24 , 55 , 29 , 0  , 66 , 68 , 63 , 9  , 26 , 31 , 64 ,
         32 , 1  , 72 , 80 , 67 , 53 , 49 , 77 , 19 , 57 , 3  , 71 , 12 , 25 ,
         42 , 7  , 11 , 41 , 60 , 28 , 93 , 8  , 76 , 15 , 91 , 82 , 46 , 61 ,
         84 , 59 , 78 , 94 , 88 , 37 , 75 , 45 , 44 , 16 , 51 , 55 , 34 , 14 ,
         62 , 50 , 38 , 73 , 73 , 65 , 43 , 90 , 22 , 81 , 69 , 92 , 23 , 74 ,
         39 , 70 , 30 , 40 , 2  , 10 , 18 , 36 , 54 , 86};


/* input */
char TitleText[MAX_TITLE_LINES][BUFF_LEN];
int TitleLines;                /* actual number of title lines */
int Natoms;                    /* number of atoms */
float Xorig,Yorig,Zorig;       /* x-, y- and z- origin */
int N1;                        /* # if incs in the slowest running direct */
float N1X1,N1Y1,N1Z1;          /* direction */
int N2;
float N2X1,N2Y1,N2Z1;          /* direction */
int N3;                        /* # if incs in the fastest running direct */
float N3X1,N3Y1,N3Z1;          /* direction */
int *IA;                       /* atomic number arraypointer */
float *Chgn;                   /* charge array pointer */
float *XC,*YC,*ZC;             /* coordinate pointers  */
float *Data;                   /* data */
int    MolecularOrbitals = 0;  /* switch to handle molecular orbitals */
int   *MolecularOrbital_p;     /* index pointer */
int    MolecularOrbitalsDefined; /* orbitals defined in file */
int    MolecularOrbital2Plot;
int    ToSave;                    /* index into the orbital array */

/* output */
float Xmax,Ymax,Zmax;
int   TypeOfData = GAUSSIAN_TYPE;
int   ProduceCoordinateFile = 1;

char InputFile[BUFF_LEN];
char OutputFile[BUFF_LEN];
char CoordinateFile[BUFF_LEN];

/* functions */
void MakeOutputFileName(char *);
int  ReadInputData(void);
int  WriteInputData(void);
int  WriteCoordinateFile(void);

/* externals */
extern char *Number2Name(int);

/**************************************************************************/
main(int args, char** argv)
/**************************************************************************/
{
    int i;

    printf("**********************************************************\n");
    printf("* Convert a 'cube' output from GaussianXX into the plot  *\n");
    printf("* format known by gOpenMol (or SCARECROW).               *\n");
    printf("*                                                        *\n");
    printf("* Leif Laaksonen (CSC) 1995                              *\n");
    printf("* Email: Leif.Laaksonen@csc.fi                           *\n");
    printf("*                               Version: 17/08/95        *\n");
    printf("**********************************************************\n\n");

    if(args < 2) {
    printf("%s",Usage);
    exit(0);}

/* go and hunt for the input switches */
    if(args > 1) {
       for(i = 1 ; i < args ; i++) {
                 /* check first that it is a '-' command */
                    if(argv[i][0] == '-' ) { /* yes it is */
/* '-m' orbital cube file                  */
      if(argv[i][1] == 'm') {
      sscanf(&argv[i][2],"%d",&MolecularOrbital2Plot);
      printf("Orbital number to plot: %d\n",MolecularOrbital2Plot);
      }
/* '-o' output file name */
      if(argv[i][1] == 'o') {
      strncpy(OutputFile,&argv[i][2],BUFF_LEN);
      }
/* '-i' input file name */
      if(argv[i][1] == 'i') {
      strncpy(InputFile,&argv[i][2],BUFF_LEN);
      }
/* '-p' prevent coordinate file output */
      if(argv[i][1] == 'p') {
       ProduceCoordinateFile = 0;
      }
     }
    }
   }

    if(InputFile[0] == '\0') {
      printf("$ERROR - supply input file name\n");
      exit(1);
    }

    printf("File names:\n");
    printf("Input file:      '%s'\n",InputFile);

    MakeOutputFileName(InputFile);
    
    printf("Output file (plot file):  '%s'\n",OutputFile);
    if(ProduceCoordinateFile)
     printf("Coordinate file (in CHARMM 'crd' format): '%s'\n",CoordinateFile);
    else
     printf("No coordinate file will be written\n");

/*  Process the data ... */

    if(ReadInputData()) {
       printf("$ERROR - can't read input data\n");
       exit(1);
    }

    if(ProduceCoordinateFile)
       (void)WriteCoordinateFile();
    (void)WriteInputData();

    printf("Job done ...\n");
}
/**************************************************************************/
void MakeOutputFileName(char *Filename)
/**************************************************************************/
{
     int i;
     int Switch;

     Switch = 1;
     if(OutputFile[0] != '\0') Switch = 0;

/* look for a dot in the name ... */
     for(i = 0 ; i < strlen(Filename) ; i++) {
         if(Filename[i] == '.') {
          if(Switch) {
           strncpy(OutputFile,Filename,(i+1));
           sprintf(&OutputFile[i+1],"plt\0");}
           strncpy(CoordinateFile,Filename,(i+1));
           sprintf(&CoordinateFile[i+1],"crd\0");
           return;
	 }
     }
       if(Switch) {
         strncpy(OutputFile,Filename,BUFF_LEN);
         strncat(OutputFile,".plt",4);}
         strncpy(CoordinateFile,Filename,BUFF_LEN);
         strncat(CoordinateFile,".crd",4);
         return;
}

/**************************************************************************/
int ReadInputData()
/**************************************************************************/
{
    FILE *Input_p;
    char  Text[BUFF_LEN];
    int   i,j,k,l,ijk,ijkl;
    int   IsInteger;
    char  Temp1[BUFF_LEN];
    char  Temp2[BUFF_LEN];
    char  Temp3[BUFF_LEN];
    char  Temp4[BUFF_LEN];
    int   Hit;
    int   tatoms;
    float txorig;
    float tyorig;
    float tzorig;

    Input_p = fopen(InputFile , "r");
    if(Input_p == NULL) {
      printf("$ERROR - can't open input file '%s'\n",InputFile);
      return(1);
    }

   printf("\nTitle in file (job title):\n");
/* first comes an unknow number of title lines (MAX lines is 5) */
    TitleLines = 0;
    for(i = 0 ; i < MAX_TITLE_LINES ; i++) {
       fgets(Text,BUFF_LEN,Input_p);
       sscanf(Text,"%s %s %s %s",Temp1,Temp2,Temp3,Temp4);
/* if Temp1 is an integer and Temp2,temp3,temp4 are floats then the title
   lines are all read */
   IsInteger = 1;
   for(j = 0 ; j < strlen(Temp1) ; j++) {
      if(isalpha(Temp1[j])) {
         IsInteger = 0;
         break;
      }
   }
   if(!IsInteger) {
       printf("%s",Text);
       strncpy(TitleText[i],Text,BUFF_LEN);
       TitleText[i][strlen(TitleText[i]) - 1] = '\0';
       TitleLines++;}
   else
       break;
   }

/* now starts the data ... */
   sscanf(Text,"%d %f %f %f",&Natoms,&Xorig,&Yorig,&Zorig);
    Xorig *= BOHR_RADIUS;
     Yorig *= BOHR_RADIUS;
      Zorig *= BOHR_RADIUS;
    printf("Number of atoms: %d, x-, y-, z-origin (in Anstrom): %f,%f,%f\n",
            abs(Natoms),Xorig,Yorig,Zorig);

   MolecularOrbitals                = 0;
   if(Natoms < 0) MolecularOrbitals = 1;

   Natoms = abs(Natoms);

   fscanf(Input_p,"%d %f %f %f",&N1,&N1X1,&N1Y1,&N1Z1);
    N1X1 *= BOHR_RADIUS;
     N1Y1 *= BOHR_RADIUS;
      N1Z1 *= BOHR_RADIUS;
    printf("Number of points: %d, in direction (x,y,z) %f %f %f\n",
            N1,N1X1,N1Y1,N1Z1);

/* check that this is a 'pure' x-coordinate */
   if(Rabs(N1X1) < SMALL) {
     printf("$ERROR - most likely your step in the x-direction (%f) is too small\n",N1X1);
     exit(20);}
   if(Rabs(N1Y1) > SMALL || Rabs(N1Z1) > SMALL) {
     printf("$ERROR - first input has to be pure x-axis (y: %f , z: %f)\n",
            N1Y1,N1Z1);
     exit(21);}

   fscanf(Input_p,"%d %f %f %f",&N2,&N2X1,&N2Y1,&N2Z1);
    N2X1 *= BOHR_RADIUS;
     N2Y1 *= BOHR_RADIUS;
      N2Z1 *= BOHR_RADIUS;
    printf("Number of points: %d, in direction (x,y,z) %f %f %f\n",
            N2,N2X1,N2Y1,N2Z1);

/* check that this is a 'pure' y-coordinate */
   if(Rabs(N2Y1) < SMALL) {
     printf("$ERROR - most likely your step in the y-direction (%f) is too small\n",N2Y1);
     exit(22);}
   if(Rabs(N2X1) > SMALL || Rabs(N2Z1) > SMALL) {
     printf("$ERROR - second input has to be pure y-axis (x: %f , z: %f)\n",
            N2X1,N2Z1);
     exit(23);}

   fscanf(Input_p,"%d %f %f %f",&N3,&N3X1,&N3Y1,&N3Z1);
    N3X1 *= BOHR_RADIUS;
     N3Y1 *= BOHR_RADIUS;
      N3Z1 *= BOHR_RADIUS;
    printf("Number of points: %d, in direction (x,y,z) %f %f %f\n",
            N3,N3X1,N3Y1,N3Z1);

/* check that this is a 'pure' z-coordinate */
   if(Rabs(N3Z1) < SMALL) {
     printf("$ERROR - most likely your step in the z-direction (%f) is too small\n",N3Z1);
     exit(24);}
   if(Rabs(N3X1) > SMALL || Rabs(N3Y1) > SMALL) {
     printf("$ERROR - first input has to be pure z-axis (x: %f , y: %f)\n",
            N3X1,N3Y1);
     exit(25);}

   IA = (int *)malloc(sizeof(int) * Natoms);
     if(IA == NULL) exit(10);
   Chgn = (float *)malloc(sizeof(float) * Natoms);
     if(Chgn == NULL) exit(11);
   XC = (float *)malloc(sizeof(float) * Natoms);
     if(XC == NULL) exit(12);
   YC = (float *)malloc(sizeof(float) * Natoms);
     if(YC == NULL) exit(13);
   ZC = (float *)malloc(sizeof(float) * Natoms);
     if(ZC == NULL) exit(14);

/* atoms ... */
   printf("Atoms...\n");
   for(i = 0 ; i < Natoms ; i++)  {
    fscanf(Input_p,"%d %f %f %f %f",&IA[i],&Chgn[i],&XC[i],&YC[i],&ZC[i]);
     XC[i] *= BOHR_RADIUS;
      YC[i] *= BOHR_RADIUS;
       ZC[i] *= BOHR_RADIUS;
    printf("Atomic number: %d, charge: %f, coord (x,y,z): %f %f %f\n",
            IA[i],Chgn[i],XC[i],YC[i],ZC[i]);

   }

   if(MolecularOrbitals) {    /* 1 */
/* molecular orbitals to handle ... */
   fscanf(Input_p,"%d",&MolecularOrbitalsDefined);

   MolecularOrbital_p = (int *)malloc(sizeof(int) * MolecularOrbitalsDefined);
   if(MolecularOrbital_p == NULL) exit(15);

   for(i = 0 ; i < MolecularOrbitalsDefined ; i++) {
       fscanf(Input_p,"%d",&MolecularOrbital_p[i]);
   }

   Hit = 0;
   for(i = 0 ; i < MolecularOrbitalsDefined ; i++) {
       if(MolecularOrbital_p[i] == MolecularOrbital2Plot) {
       Hit = MolecularOrbital2Plot;
       ToSave = i;
       break;}
   }

   if(!Hit) {
   printf("$ERROR - Specified orbital nr: %d is not in the file\n",
          MolecularOrbital2Plot);
   exit(16);}

   Data = (float *)malloc(sizeof(float) * N1 * N2 * N3 * 
                                          MolecularOrbitalsDefined);
     if(Data == NULL) exit(17);

   for(i = 0   ; i < N1 ; i++)  { /* 2 */
    for(j = 0  ; j < N2 ; j++)  { /* 3 */
     for(k = 0 ; k < N3 ; k++)  { /* 4 */
      for(l = 0; l < MolecularOrbitalsDefined ; l++) { /* 5 */

     ijkl = i + N1 * j + N1 * N2 * k + N1 * N2 * N3 * l;

     ijk = fscanf(Input_p,"%f",&Data[ijkl]);

     if(ijk != 1) {
       printf("$ERROR - in reading the grid data\n");
       printf("$ERROR - at ijkl: %d %d %d %d\n",i,j,k,l);
       exit(18);
     }
       }/* 5 */
      } /* 4 */
     }  /* 3 */
    }   /* 2 */
   } /* end *1* */
   else {  /* start *1* */

   Data = (float *)malloc(sizeof(float) * N1 * N2 * N3);
     if(Data == NULL) exit(19);

   for(i = 0   ; i < N1 ; i++)  { /* 2 */
    for(j = 0  ; j < N2 ; j++)  { /* 3 */
     for(k = 0 ; k < N3 ; k++)  { /* 4 */

     ijk = i + N1 * j + N1 * N2 * k;

     ijkl = fscanf(Input_p,"%f",&Data[ijk]);
     if(ijkl != 1) {
       printf("$ERROR - in reading the grid data\n");
       printf("$ERROR - at ijk: %d %d %d\n",i,j,k);
       exit(20);
     }
     } /* 4 */
    }  /* 3 */
   }   /* 2 */
  } /* end *1* */
   fclose(Input_p);

   return(0);
}

/**************************************************************************/
int WriteInputData()
/**************************************************************************/
{
    FILE *Output_p;
    int   i,j,k,l,ijk,ijkl;
    int Items;
    float Help1;
    float Help2;

    Output_p = fopen(OutputFile,"w");
    if(Output_p == NULL) {
      printf("$ERROR - can't open output file: '%s' \n",OutputFile);
      return(1);
    }

    i = 3;
    FWRITE(&i,sizeof(int));
    FWRITE(&TypeOfData , sizeof(int));

      FWRITE(&N3   , sizeof(int));
       FWRITE(&N2  , sizeof(int));
        FWRITE(&N1 , sizeof(int));

      Help1  = Zorig;
       Help2 = Zorig + (N3 - 1) * N3Z1;
        FWRITE(&Help1 , sizeof(float));
        FWRITE(&Help2 , sizeof(float));
      Help1  = Yorig;
       Help2 = Yorig + (N2 - 1) * N2Y1;
        FWRITE(&Help1 , sizeof(float));
        FWRITE(&Help2 , sizeof(float));
      Help1  = Xorig;
       Help2 = Xorig + (N1 - 1) * N1X1;
        FWRITE(&Help1 , sizeof(float));
        FWRITE(&Help2 , sizeof(float));

   if(MolecularOrbitals) {         

   for(k = 0   ; k < N3 ; k++)  {  /* 2 */
    for(j = 0  ; j < N2 ; j++)  {  /* 3 */
     for(i = 0 ; i < N1 ; i++)  {  /* 4 */

     ijk = i + N1 * j + N1 * N2 * k + N1 * N2 * N3 * ToSave;

     FWRITE(&Data[ijk] , sizeof(float));
     }                             /* 4 */
    }                              /* 3 */
   }                               /* 2 */
   }
   else {
   for(k = 0   ; k < N3 ; k++)  {
    for(j = 0  ; j < N2 ; j++)  { 
     for(i = 0 ; i < N1 ; i++)  {

     ijk = i + N1 * j + N1 * N2 * k;

     FWRITE(&Data[ijk] , sizeof(float));
     }
    }
   }
 }
   fclose(Output_p);

   return(0);
}

/**************************************************************************/
int WriteCoordinateFile()
/**************************************************************************/
{
    int   i,j;
    FILE *coord_p;
    char  AtomName[2];

    coord_p = fopen(CoordinateFile,"w");
    if(coord_p == NULL) {
      printf("$ERROR - can't open output file: '%s'",CoordinateFile);
      exit(1);
    }

    fprintf(coord_p,"* ++Automatic output generated from Gaussian++\n");
    fprintf(coord_p,"* ++using the 'cube' keyword                ++\n");

    for(i = 0 ; i < TitleLines ; i++) 
        fprintf(coord_p,"* %.70s\n",TitleText[i]);

    fprintf(coord_p,"*  \n");
    fprintf(coord_p,"%5d \n",Natoms);

    for(i = 0 ; i < Natoms ; i++) {

      for(j = 0 ; j < strlen(AtomSymbols)/4 ; j++) {
         if(IA[i] == AtomSymbol_p[j]) { 
         strncpy(AtomName,&AtomSymbols[4 * j],2);
         AtomName[2] = '\0';
         break;}
      }
      fprintf(coord_p,
      "%5d%5d %-4.4s %-4.4s%10.5f%10.5f%10.5f %4.4s %-4d%10.5f \n",
      (i+1),1,"GAUS",AtomName,XC[i],YC[i],ZC[i],"GAUS",1,0.0);
      }

    fclose(coord_p);

    return(0);

}