/*******************************************************************************
*                                                                              *
*                                   Viewmol                                    *
*                                                                              *
*                                I N P U T . C                                 *
*                                                                              *
*                 Copyright (c) Joerg-R. Hill, October 2003                    *
*                                                                              *
********************************************************************************
*
* $Id: input.c,v 1.7 2004/08/29 14:52:57 jrh Exp $
* $Log: input.c,v $
* Revision 1.7  2004/08/29 14:52:57  jrh
* Release 2.4.1
*
* Revision 1.6  2003/11/07 11:05:02  jrh
* Release 2.4
*
* Revision 1.5  2000/12/10 15:09:24  jrh
* Release 2.3
*
* Revision 1.4  1999/05/24 01:25:54  jrh
* Release 2.2.1
*
* Revision 1.3  1999/02/07 21:50:31  jrh
* Release 2.2
*
* Revision 1.2  1998/01/26 00:48:07  jrh
* Release 2.1
*
* Revision 1.1  1996/12/10  18:41:26  jrh
* Initial revision
*
*/
#include<ctype.h>
#include<locale.h>
#include<stdlib.h>
#include<stdio.h>
#include<string.h>
#include<unistd.h>
#include<X11/Intrinsic.h>
#include<Xm/Xm.h>
#include "viewmol.h"
#include "dialog.h"

#define MAX(a, b) (a) > (b) ? (a) : (b)
#define MIN(a, b) (a) < (b) ? (a) : (b)

int input(int, char **);
int load(Widget, caddr_t, XmAnyCallbackStruct *);
char *findFileType(char *, char *);
int  readFile(char *, char *);
double *readMOs(struct MOLECULE *, FILE *, int, int *, char *);
int  readOccupation(struct MOLECULE *, FILE *, int, char *);
void relint(struct MOLECULE *);
int  assignBasisset(struct MOLECULE *, int, char *);
int  moveBasissetData(struct MOLECULE *, struct BASISSET *, int, int);
void changeExponentPointers(struct BASISSET *, long, int);
void sortBasis(struct MOLECULE *);
int  compar(const void *, const void *);
struct SAVE *saveGeometry(struct MOLECULE *);
void wrbas(struct MOLECULE *);
struct MOLECULE *initMolecule(int);
void deleteMolecule(Widget, caddr_t, XmAnyCallbackStruct *);
void reformat(char *, char *);
PyMoleculeSpecObject *molecule_new(void);

extern void GetMessageBoxButton(Widget, XtPointer, caddr_t);
extern char *selectFile(char *, char *, int);
extern char *getStringResource(Widget, char *);
extern int messgb(Widget, int, char *, struct PushButtonRow *, int);
extern void *getmem(size_t, size_t);
extern void *expmem(void *, size_t, size_t);
extern void fremem(void **);
extern int makeConnectivity(struct MOLECULE *, int, int);
extern void makeUnitCell(double, double, double, double, double, double, double *, int);
extern void shiftUnitCell(struct MOLECULE *);
extern int getNumber(char *, int *);
extern void cnorm(int);
extern void modcmo(struct MOLECULE *, int);
extern void assignBasisName(int);
extern int checkFile(char **);
extern void setWindowTitle(Widget, char *);
extern void setAtom(struct MOLECULE *);
extern void setMenuItem(int, int);
extern void setMenu(struct MOLECULE *);
extern void redraw(int);
extern void inert(struct MOLECULE *);
extern void setBond(struct MOLECULE *);
extern Widget makeViewerMenu(Widget);
extern void makeSpectrumMenu(void);
extern void makeMODiagramMenu(void);
extern void makeHistoryMenu(void);
extern void setWindow(int);
extern void showTitle(Widget,  caddr_t, XmDrawingAreaCallbackStruct *);
extern void drawMolecule(Widget, caddr_t, XmDrawingAreaCallbackStruct *);
extern void quitSpectrum(Widget, caddr_t, XmAnyCallbackStruct *);
extern void quitHistory(Widget, caddr_t, XmAnyCallbackStruct *);
extern void quitMODiagram(Widget, caddr_t, XmAnyCallbackStruct *);
extern void selectMolecule(Widget, XtPointer, XmToggleButtonCallbackStruct *);
extern void runScript(Widget, caddr_t, XmAnyCallbackStruct *);
extern void readGrid(struct GRIDOBJECT *, FILE *, char *);
extern void adjustMONumbers(struct MOLECULE *);
extern void addBoundaryAtoms(struct MOLECULE *);
extern void getSumFormula(struct MOLECULE *, char *, int);

extern struct MOLECULE *molecules;
extern struct OPTION *options;
extern struct WINDOW windows[];
extern float *transObject, *rotObject;
extern int nopt, nmolecule, imol, moveItem;
extern int showWarning, debug;
extern Widget topShell;

int input(int argc, char **argv)
{
  static struct PushButtonRow buttons1[] = {{"exit", GetMessageBoxButton,
                                             (XtPointer)0, NULL}};
  const char *locale;
  char command[MAXLENLINE], line[MAXLENLINE];
  char *word=NULL;
  int readNext, n;
  register int i;

  if (argc > 1)
  {
    /* we got an option */
    n=1;
    while (n < argc)
    {
      if (*argv[n] == '-')
      {
        argv[n]++;
        readNext=TRUE;
        for (i=0; i<nopt; i++)
        {
          if (!strncmp(options[i].flag, argv[n], 8))
          {
            sprintf(command, options[i].command, argv[n+1]);
            readNext=FALSE;
            break;
          }
        }
        if (readNext)
        {
          word=getStringResource(topShell, "unknownParameter");
          sprintf(line, word, --argv[n]);
          messgb(topShell, 3, line, buttons1, 1);
          exit(-1);
        }
        word=argv[n+1];
        locale=setlocale(LC_ALL, "C");
        i=readFile(command, word); 
        (void)setlocale(LC_ALL, locale);
        if (!i) return(FALSE);
        n+=2;
      }
      else
      {
        /* no option, but hopefully a file name */
        if ((word=findFileType(command, argv[n])) == NULL) exit(-1);
        if (!strncmp(command, "python", 6))
          runScript((Widget)0, (caddr_t)argv[n], (XmAnyCallbackStruct *)0);
        else
        {
          locale=setlocale(LC_ALL, "C");
          i=readFile(command, word);
          if (strcmp(word, argv[n])) unlink(word);
          (void)setlocale(LC_ALL, locale);
          if (!i) return(FALSE);
        }
        n++;
      }
    }
  }
  else  /* argv <= 1 */
  {
    readNext=TRUE;
    for (i=0; i<nopt; i++)
    {
      if (!strcmp(options[i].flag, "default"))
      {
        strcpy(command, options[i].command);
        readNext=FALSE;
        break;
      }
    }
    if (readNext) return(FALSE);
    locale=setlocale(LC_ALL, "C");
    i=readFile(command, word);
    (void)setlocale(LC_ALL, locale);
    if (!i) return(FALSE);
  }
  if (nmolecule == 0)
    return(FALSE);
  else
    return(TRUE);
}

int load(Widget w, caddr_t data, XmAnyCallbackStruct *dummy)
{
  XmToggleButtonCallbackStruct button;
  static char *filename="";
  const char *locale;
  char line[MAXLENLINE], *f;
  int i;
  register int ret;

  if ((char *)data == NULL)
    filename=selectFile("*", filename, FALSE);
  else
    filename=(char *)data;
  if (filename != NULL)
  {
    if ((f=findFileType(line, filename)) == NULL) return(FALSE);
    windows[VIEWER].set=nmolecule;
    locale=setlocale(LC_ALL, "C");
    ret=readFile(line, filename);
    if (strcmp(f, filename)) unlink(f);
    (void)setlocale(LC_ALL, locale);
    if (!ret) return(FALSE);
    windows[VIEWER].set=imol;
    setAtom(&molecules[imol]);
    if (!molecules[imol].title[0]) getSumFormula(&molecules[imol], molecules[imol].title, FALSE);
    /* Check uniqueness of title */
    for (i=0; i<nmolecule; i++)
    {
      if (i != imol)
      {
        if (!strcmp(molecules[i].title, molecules[imol].title))
        {
          sprintf(line, "%s <%d>", molecules[imol].title, nmolecule);
          strncpy(molecules[imol].title, line, MAXLENLINE);
        break;
        }
      }
    }
    if (molecules[imol].unitcell) addBoundaryAtoms(&molecules[imol]);
    makeConnectivity(&molecules[imol], TRUE, TRUE);
    inert(&molecules[imol]);
    moveItem=imol+MOLECULES;
    setMenuItem(VIEWER1_SAVE, True);
    setMenuItem(VIEWER1_DELETE, True);
    setMenuItem(VIEWER1_BUILD, True);
    windows[VIEWER].menu=makeViewerMenu(windows[VIEWER].widget);
    setWindow(imol);
    button.set=TRUE;
    if (windows[SPECTRUM].widget != NULL && nmolecule > 1)
    {
      if (molecules[imol].nmodes > 0)
      {
        windows[SPECTRUM].set=imol;
        makeSpectrumMenu();
        (void)selectMolecule((Widget)0, (XtPointer)&windows[SPECTRUM].selectMenu[imol],
                             &button);
      }
    }
    if (windows[MO].widget != NULL && nmolecule > 1)
    {
      if (molecules[imol].nbasfu > 0)
      {
        windows[MO].set=imol;
        makeMODiagramMenu();
        (void)selectMolecule((Widget)0, (XtPointer)&windows[MO].selectMenu[imol],
                             &button);
      }
    }
    if (windows[HISTORY].widget != NULL && nmolecule > 1)
    {
      if (molecules[imol].nhist > 0)
      {
        windows[HISTORY].set=imol;
        makeHistoryMenu();
        (void)selectMolecule((Widget)0, (XtPointer)&windows[HISTORY].selectMenu[imol],
                             &button);
      }
    }
    XtRemoveAllCallbacks(windows[VIEWER].widget, XmNexposeCallback);
    XtAddCallback(windows[VIEWER].widget, XmNexposeCallback, (XtCallbackProc)drawMolecule, NULL);
    redraw(VIEWER);
    return(TRUE);
  }
  return(FALSE);
}

char *findFileType(char *command, char *filename)
{
  FILE *file;
  static struct PushButtonRow buttons[] = {{"continue", GetMessageBoxButton, (XtPointer)0, NULL}};
  static char *temp;
  char *buffer, *word;
  register int i;

  if ((file=fopen(filename, "r")) == NULL)
  {
    word=getStringResource(topShell, "noFile");
    sprintf(command, word, filename);
    messgb(topShell, 3, command, buttons, 1);
    return(NULL);
  }
  buffer=(char *)getmem((size_t)1025, sizeof(char));
  fread((void *)buffer, sizeof(char), (size_t)1024, file);
  fclose(file);
  buffer[1024]='\0';
  for (i=0; i<nopt; i++)
  {
    if (options[i].identifier && strstr(buffer, options[i].identifier))
    {
      if (strstr(buffer, "\r\n"))  /* we got a DOS formatted file */
      {
        temp=tmpnam(NULL);
        reformat(filename, temp);
        sprintf(command, options[i].command, temp);
        return(temp);
      }
      else
      {
        sprintf(command, options[i].command, filename);
        return(filename);
      }
      /* NOTREACHED */
      break;
    }
  }
  fremem((void **)&buffer);
  word=getStringResource(topShell, "formatNotRecognized");
  sprintf(command, word, filename);
  messgb(topShell, 3, command, buttons, 1);
  return(NULL);
}

int readFile(char *command, char *filename)
{
  static struct PushButtonRow buttons1[] = {{"exit", GetMessageBoxButton,
                                             (XtPointer)0, NULL}};
  static struct PushButtonRow buttons2[] = {{"exit", GetMessageBoxButton,
                                             (XtPointer)0, NULL},
                                            {"continue", GetMessageBoxButton,
                                             (XtPointer)1, NULL}};
  FILE *file;
  struct MOLECULE *mol;
  struct BASISSET *oldbasisset;
  char line[MAXLENLINE], label[MAXLENLINE];
  char fixed[4], *word, *p, *w;
  GLdouble x, y, z;
  double *cmo=NULL, *oldexponents;
  double a=0.0, b=0.0, c=0.0, alpha=0.0, beta=0.0, gamma=0.0;
  double aopt=0.0, bopt=0.0, copt=0.0, alphaopt=0.0, betaopt=0.0, gammaopt=0.0;
  double radfac=-1.0, radfac1;
  double bravaisMatrix[9];
  size_t mna=20, mnn=60, mnh=50, mns, mnp;
  static int shell=0, primitive=0;
  int readNext=FALSE, ncoeff=0, *help;
  int frac=FALSE, ierror=FALSE;
  int namax=0, istar;
  register int i, j, k, m;

  word=command;
  if (!checkFile(&word))
  {
    word=getStringResource(topShell, "noFile");
    strtok(command, " \t");
    sprintf(line, word, command);
    messgb(topShell, 3, line, buttons1, 1);
    return(FALSE);
  }
  if ((file=popen(command, "r")) == NULL)
  {
    word=getStringResource(topShell, "cannotExecute");
    sprintf(line, word, command);
    messgb(topShell, 3, line, buttons1, 1);
    return(FALSE);
  }

  mol=initMolecule(mna);
  word=fgets(line, MAXLENLINE, file);
  while (word != NULL)
  {
    if (strstr(line, "$coord"))
    {
      i=0;
      word=strtok(line, " \t");
      if ((word=strtok(NULL, " \t")) != NULL)
      {
        if (!strcmp(word, "fractional"))
        {
          frac=TRUE;
          if (radfac < 0.0)
          {
            if ((word=strtok(NULL, " \t")) != NULL)
              radfac=atof(word);
            else
              radfac=1.0;
          }
        }
        else
          if (radfac < 0.0) radfac=atof(word);
      }
      else
        if (radfac < 0.0) radfac=1.0;
      while (fgets(line, MAXLENLINE, file) != NULL)
      {
        if (strchr(line, '$')) break;
        strcpy(fixed, "   ");
        sscanf(line, "%lf%lf%lf%s%s", &x, &y, &z, mol->atoms[i].name, fixed);
        mol->atoms[i].x=x*radfac;
        mol->atoms[i].y=y*radfac;
        mol->atoms[i].z=z*radfac;
        mol->atoms[i].name[0]=toupper(mol->atoms[i].name[0]);
        mol->atoms[i].name[1]=tolower(mol->atoms[i].name[1]);
        mol->atoms[i].ref=i;
        mol->atoms[i].flags=ORIGINAL;
        if (strpbrk(fixed, "xX")) mol->atoms[i].flags|=X_FIXED;
        if (strpbrk(fixed, "yY")) mol->atoms[i].flags|=Y_FIXED;
        if (strpbrk(fixed, "zZ")) mol->atoms[i].flags|=Z_FIXED;
	mol->atoms[i].basis=NULL;
	mol->atoms[i].basisname[0]='\0';
        if (++i >= mna)
        {
          mna+=20;
          mol->atoms=(struct ATOM *)expmem((void *)mol->atoms, mna,
                                           sizeof(struct ATOM));
        }
      }
      mol->na=i;
      mna=i;
      mol->atoms=(struct ATOM *)expmem((void *)mol->atoms, mol->na,
                                       sizeof(struct ATOM));
      readNext=FALSE;
    }
    else if (strstr(line, "$title"))
    {
      fgets(line, MAXLENLINE, file);
      if (line[0] != '$')
      {
        word=line+strlen(line)-1;
        while (*word == '\n' || *word == ' ')
        {
          *word='\0';
          word--;
        }
        if (line[0] != '\0') strcpy(mol->title, line);
        readNext=TRUE;
      }
      else
        readNext=FALSE;
    }
    else if (strstr(line, "$vibrational normal modes"))
    {
      mol->cnm=(double *)getmem(mnn*mnn, sizeof(double));
      i=0;
      while (fgets(line, MAXLENLINE, file) != NULL)
      {
        if (strchr(line, '$')) break;
        word=strtok(line, " \t");
	if (strlen(word) < 3)
          word=strtok(NULL, " \t");
        while ((word=strtok(NULL, " \t")) != NULL)
        {
          mol->cnm[i]=atof(word); 
          if (++i >= mnn*mnn)
          {
            mnn+=60;
            mol->cnm=(double *)expmem((void *)mol->cnm, mnn*mnn, sizeof(double));
          }
        }
      }
      readNext=FALSE;
    }
    else if (strstr(line, "$vibrational spectrum"))
    {
      mol->normal_modes=(struct NORMAL_MODE *)getmem(mnn,
                                                sizeof(struct NORMAL_MODE));
      i=0;
      while (fgets(line, MAXLENLINE, file) != NULL)
      {
        if (strchr(line, '$')) break;
        sscanf(line, "%s%lf%lf%lf", mol->normal_modes[i].sym,
                                    &(mol->normal_modes[i].wavenumber),
                                    &(mol->normal_modes[i].ir_intensity),
                                    &(mol->normal_modes[i].raman_intensity));
        if (++i >= mnn)
        {
          mnn+=60;
          mol->normal_modes=(struct NORMAL_MODE *)expmem((void *)mol->normal_modes,
                                                  mnn, sizeof(struct NORMAL_MODE));
        }
      }
      mol->nmodes=i;
      if (mol->nmodes > 0)
      {
        mol->normal_modes=(struct NORMAL_MODE *)expmem((void *)mol->normal_modes,
                           mol->nmodes, sizeof(struct NORMAL_MODE));
        relint(mol);
      }
      readNext=FALSE;
    }
    else if (strstr(line, "$grad") && !strstr(line, "$grad_"))
    {
      if (radfac < 0.0)
      {
        (void)strtok(line, " \t\n");
        if ((word=strtok(NULL, " \t\n")) != NULL)
        {
          radfac=strtod(word, &p);
          if (p == word) radfac=1.0;
        }
        else
          radfac=1.0;
      }
      mol->history=(struct COORDS *)getmem(mnh*mna, sizeof(struct COORDS));
      mol->optimization=(struct OPTIMIZATION *)getmem(mnh, sizeof(struct OPTIMIZATION));
      i=j=k=0;
      namax=0;
      m=0;
      while (fgets(line, MAXLENLINE, file) != NULL)
      {
        if (strchr(line, '$')) break;
        if (strstr(line, "cycle"))
        {
          strtok(line, "=");
          strtok(NULL, "=");
          mol->optimization[i].energy=atof(strtok(NULL, "="));
          mol->optimization[i].gnorm=atof(strtok(NULL, "="));
          mol->optimization[i].coords=j;
          if (namax < m)
          {
            namax=m;
            mol->cycle=i;
          }
          if (i > 0) mol->optimization[i-1].natoms=m;
          m=0;
          if (++i >= mnh)
          {
            mnh+=50;
            mol->history=(struct COORDS *)expmem((void *)mol->history,
                                         mnh*mna, sizeof(struct COORDS));
            mol->optimization=(struct OPTIMIZATION *)expmem((void *)mol->optimization,
                                              mnh, sizeof(struct OPTIMIZATION));
          }
        }
        else if (strstr(line, "unitcell"))
        {
          if (strstr(line, "vectors"))
          {
            fgets(line, MAXLENLINE, file);
            bravaisMatrix[0]=radfac*atof(strtok(line, " \t"));
            bravaisMatrix[1]=radfac*atof(strtok(NULL, " \t"));
            bravaisMatrix[2]=radfac*atof(strtok(NULL, " \t"));
            fgets(line, MAXLENLINE, file);
            bravaisMatrix[3]=radfac*atof(strtok(line, " \t"));
            bravaisMatrix[4]=radfac*atof(strtok(NULL, " \t"));
            bravaisMatrix[5]=radfac*atof(strtok(NULL, " \t"));
            fgets(line, MAXLENLINE, file);
            bravaisMatrix[6]=radfac*atof(strtok(line, " \t"));
            bravaisMatrix[7]=radfac*atof(strtok(NULL, " \t"));
            bravaisMatrix[8]=radfac*atof(strtok(NULL, " \t"));
            aopt=-1.0;
          }
          else
          {
            word=strtok(line, " \t");
            aopt=atof(strtok(NULL, " \t"));
            bopt=atof(strtok(NULL, " \t"));
            copt=atof(strtok(NULL, " \t"));
            alphaopt=atof(strtok(NULL, " \t"));
            betaopt=atof(strtok(NULL, " \t"));
            gammaopt=atof(strtok(NULL, " \t"));
            mol->optimization[i-1].coords=j;
          }
          makeUnitCell(radfac*aopt, radfac*bopt, radfac*copt, alphaopt, betaopt,
                       gammaopt, bravaisMatrix, FALSE);
        }
        else
        {
          sscanf(line, "%lf%lf%lf", &x, &y, &z);
          strtok(line, " \t");
          strtok(NULL, " \t");
          strtok(NULL, " \t");
          if (strtok(NULL, " \t") != NULL)
          {
            mol->history[j].x=x;
            mol->history[j].y=y;
            mol->history[j].z=z;
            m++;
            if (++j >= mnh*mna)
            {
              mna+=20;
              mol->history=(struct COORDS *)expmem((void *)mol->history,
                                           mnh*mna, sizeof(struct COORDS));
            }
          }
          else
          {
            mol->history[k].gx=x;
            mol->history[k].gy=y;
            mol->history[k++].gz=z;
          }
        }
      }
      if (i > 0)
      {
        mol->nhist=i;
        mol->cycle=i;
        mol->optimization[i-1].natoms=m;
        if (namax < m)
        {
          namax=m;
          mol->cycle=i;
        }
        mol->atoms=(struct ATOM *)expmem((void *)mol->atoms, namax+1,
                                         sizeof(struct ATOM));
        if (j > 0)
          mol->history=(struct COORDS *)expmem((void *)mol->history, j,
                                               sizeof(struct COORDS));
        if (mol->nhist > 0)
          mol->optimization=(struct OPTIMIZATION *)expmem((void *)mol->optimization,
                                            mol->nhist, sizeof(struct OPTIMIZATION));
        mol->emax=mol->emin=mol->optimization[0].energy;
        mol->gmax=mol->optimization[0].gnorm;
        for (i=1; i<mol->nhist; i++)
        {
          mol->emax=MAX(mol->emax, mol->optimization[i].energy);
          mol->emin=MIN(mol->emin, mol->optimization[i].energy);
          mol->gmax=MAX(mol->gmax, mol->optimization[i].gnorm);
        }
        if (mol->gmax == 0.0) windows[HISTORY].mode=SCALES | ENERGY;
      }
      else
      {
        fremem((void **)&mol->history);
        fremem((void **)&mol->optimization);
      }
      readNext=FALSE;
    }
    else if (strstr(line, "$scfmo"))
    {
      if (strstr(line, "symmetrized")) mol->needMoloch=TRUE;
      if (strstr(line, "gaussian")) mol->gaussian=TRUE;
      cmo=readMOs(mol, file, ALPHAANDBETA, &ncoeff, line);
      readNext=FALSE;
    }
    else if (strstr(line, "$uhfmo_alpha"))
    {
      if (strstr(line, "symmetrized")) mol->needMoloch=TRUE;
      if (strstr(line, "gaussian")) mol->gaussian=TRUE;
      cmo=readMOs(mol, file, ALPHA, &ncoeff, line);
      readNext=FALSE;
    }
    else if (strstr(line, "$uhfmo_beta"))
    {
      if (strstr(line, "symmetrized")) mol->needMoloch=TRUE;
      if (strstr(line, "gaussian")) mol->gaussian=TRUE;
      cmo=readMOs(mol, file, BETA, &ncoeff, line);
      readNext=FALSE;
    }
    else if (strstr(line, "$grid"))
    {
      if (mol->gridObjects == NULL)
        mol->gridObjects=(struct GRIDOBJECT *)getmem((size_t)1, sizeof(struct GRIDOBJECT));
      else
        mol->gridObjects=(struct GRIDOBJECT *)expmem((void *)mol->gridObjects, 
                         (size_t)(mol->ngridobjects+1), sizeof(struct GRIDOBJECT));
      readGrid(&mol->gridObjects[mol->ngridobjects], file, line);
      mol->ngridobjects++;
      readNext=FALSE;
    }
    else if (strstr(line, "$closed shells"))
    {
      readNext=FALSE;
      if (!readOccupation(mol, file, ALPHAANDBETA, line)) continue;
    }
    else if (strstr(line, "$alpha shells"))
    {
      readNext=FALSE;
      if (!readOccupation(mol, file, ALPHA, line)) continue;
    }
    else if (strstr(line, "$beta shells"))
    {
      readNext=FALSE;
      if (!readOccupation(mol, file, BETA, line)) continue;
    }
    else if (strstr(line, "$atoms"))
    {
      while (fgets(line, MAXLENLINE, file) != NULL)
      {
        if (strchr(line, '$')) break;
        w=strrchr(line, '\\');
        word=strtok(line, " \t");
        word=strtok(NULL, " \t");
        help=(int *)getmem(mna, sizeof(int));
        i=getNumber(word, help);
        while (w != NULL)
        {
          fgets(line, MAXLENLINE, file);
          w=strrchr(line, '\\');
          if (strstr(line, "basis") != NULL)
          {
            word=strtok(line, "=");
            word=strtok(NULL, "\\\n");
            p=word+strlen(word);
            while (*(--p) == ' ');
            k=(int)(p-word)+1;
            k=k < MAXLENBASISNAME-1 ? k : MAXLENBASISNAME-1;
            for (j=0; j<i; j++)
            {
              strncpy(mol->atoms[help[j]-1].basisname, word, k);
              mol->atoms[help[j]-1].basisname[k]='\0';
            }
          }
        }
        fremem((void **)&help);
      }
      readNext=FALSE;
    }
    else if (strstr(line, "$basis"))
    {
      ierror=FALSE;
      mns=100;
      mnp=200;
      if (mol->basisset == NULL)
        mol->basisset=(struct BASISSET *)getmem(mns, sizeof(struct BASISSET));
      if (mol->exponents == NULL)
        mol->exponents=(double *)getmem(mnp, 2*sizeof(double));
      shell=0;
      primitive=0;

      while (fgets(line, MAXLENLINE, file) != NULL)
      {
        if (strchr(line, '*'))
        {
	  newSet:
          if (shell != 0) mol->basisset[shell-1].next=(struct BASISSET *)NULL;
          fgets(line, MAXLENLINE, file);
          if (strchr(line, '$')) break;
          if (!assignBasisset(mol, shell, line))
	  {
	    /* This basis set is not assigned to any atom, skip it */
	    istar=0;
	    while (fgets(line, MAXLENLINE, file) != NULL)
	    {
	      if (strchr(line, '$')) break;
	      if (strchr(line, '*')) istar++;
	      if (istar > 1) goto newSet;
	    }
	    break;
	  }
          fgets(line, MAXLENLINE, file);
          fgets(line, MAXLENLINE, file);
        }
        mol->basisset[shell].first=&mol->exponents[primitive];
        word=strtok(line, " \t");
        mol->basisset[shell].nprim=atoi(word);
        word=strtok(NULL, " \t");
        if (isdigit(*word))
          mol->basisset[shell].ang=atoi(word);
        else
        {
          switch (*word)
          {
            case 's': mol->basisset[shell].ang=0;
                      break;
            case 'p': mol->basisset[shell].ang=1;
                      break;
            case 'd': mol->basisset[shell].ang=2;
                      break;
            default:  mol->basisset[shell].ang=(*word)-'c';
                      break;
          }
        }
        for (j=0; j<mol->basisset[shell].nprim; j++)
        {
          fgets(line, MAXLENLINE, file);
          word=strtok(line, " \t");
          mol->exponents[primitive++]=atof(word);
          word=strtok(NULL, " \t");
          mol->exponents[primitive++]=atof(word);
          if (primitive >= mnp)
          {
            mnp+=100;
            oldexponents=mol->exponents;
            mol->exponents=(double *)expmem((void *)mol->exponents, mnp, 2*sizeof(double));
            if (mol->exponents != oldexponents)
              changeExponentPointers(mol->basisset, mol->exponents-oldexponents, shell);
          }
        }
        if (++shell >= mns)
        {
          mns+=100;
          oldbasisset=mol->basisset;
          mol->basisset=(struct BASISSET *)getmem(mns, sizeof(struct BASISSET));
          ierror|=moveBasissetData(mol, oldbasisset, shell, FALSE);
          fremem((void **)&oldbasisset);
        }
        mol->basisset[shell-1].next=&mol->basisset[shell];
      } /* end of "while (fgets(line, MAXLENLINE, file) != NULL)" */
      readNext=FALSE;
      oldexponents=mol->exponents;
      if (primitive > 0)
        mol->exponents=(double *)expmem((void *)mol->exponents, primitive, sizeof(double));
      if (mol->exponents != oldexponents)
        changeExponentPointers(mol->basisset, mol->exponents-oldexponents, shell);
      oldbasisset=mol->basisset;
      if (shell > 0)
      {
        mol->basisset=(struct BASISSET *)getmem(shell, sizeof(struct BASISSET));
        ierror|=moveBasissetData(mol, oldbasisset, shell, TRUE);
        fremem((void **)&oldbasisset);
      }
      if (ierror)
      {
        fremem((void **)&mol->basisset);
        fremem((void **)&mol->exponents);
      }
      else
      {
        if (debug == 2) wrbas(mol);
        cnorm(nmolecule);
      }
    }
    else if (strstr(line, "$unitcell"))
    {
      if (strstr(line, "vectors"))
      {
        fgets(line, MAXLENLINE, file);
        bravaisMatrix[0]=atof(strtok(line, " \t"));
        bravaisMatrix[1]=atof(strtok(NULL, " \t"));
        bravaisMatrix[2]=atof(strtok(NULL, " \t"));
        fgets(line, MAXLENLINE, file);
        bravaisMatrix[3]=atof(strtok(line, " \t"));
        bravaisMatrix[4]=atof(strtok(NULL, " \t"));
        bravaisMatrix[5]=atof(strtok(NULL, " \t"));
        fgets(line, MAXLENLINE, file);
        bravaisMatrix[6]=atof(strtok(line, " \t"));
        bravaisMatrix[7]=atof(strtok(NULL, " \t"));
        bravaisMatrix[8]=atof(strtok(NULL, " \t"));
        a=-1.0;
      }
      else
      {
        word=strtok(line, " \t");
        a=atof(strtok(NULL, " \t"));
        b=atof(strtok(NULL, " \t"));
        c=atof(strtok(NULL, " \t"));
        alpha=atof(strtok(NULL, " \t"));
        beta=atof(strtok(NULL, " \t"));
        gamma=atof(strtok(NULL, " \t"));
      }
      readNext=TRUE;
    }
    else if (strstr(line, "$pople"))
    {
      if (strstr(line, "6d"))  mol->pured=FALSE;
      if (strstr(line, "10f")) mol->puref=FALSE;
      if (strstr(line, "15g")) mol->pureg=FALSE;
      readNext=TRUE;
    }
    else if (strstr(line, "$symmetry"))
    {
      word=strtok(line, " \t\n");
      if ((word=strtok(NULL, " \t\n")) != NULL)
      {
        strcpy(mol->pgroup, word);
        word=mol->pgroup;
        while (*word && word < mol->pgroup+5)
        {
          *word=tolower(*word);
          word++;
        }
      }
      readNext=TRUE;
    }
    else if (strstr(line, "$error"))
    {
      word=strtok(line, " \t");
      word=strtok(NULL, " \t");
      if ((w=getStringResource(topShell, word)) == NULL)
      {
        w=getStringResource(topShell, "unknownErrorMessage");
        i=TRUE;
      }
      else
      {
	/* Get severity of error */
        i=atoi(strtok(NULL, " \t"));
        word=strtok(NULL, " \t\n");
      }
      sprintf(label, w, word);
      if (i)
      {
	if (showWarning) messgb(topShell, 3, label, buttons1, 1);
        return(FALSE);
      }
      else
      {
        if (showWarning)
          if (!messgb(topShell, 3, label, buttons2, 2)) return(FALSE);
      }
      readNext=TRUE;
    }
    else if (strstr(line, "$end"))
    {
      pclose(file);
      if (mol->basisset != NULL)
      {
        if (!mol->gaussian) sortBasis(mol);
        assignBasisName(nmolecule);
      }
      if (mol->orbitals != NULL && ncoeff != 0)
      {
        if (mol->basisset != NULL && !mol->needMoloch)
        {
          modcmo(mol, mol->gaussian);
          fremem((void **)&cmo);
        }
      }
      if (mol->orbitals != NULL && mol->ngridobjects > 0)
        adjustMONumbers(mol);
      if (a != 0.0)
      {
        if (a < 0.0)
        {
          for (i=0; i<9; i++)
            bravaisMatrix[i]*=radfac;
        }
        makeUnitCell(a*radfac, b*radfac, c*radfac, alpha, beta, gamma, bravaisMatrix, frac);
      for (i=0; i<mol->unitcell->nc; i++)
        {
          mol->unitcell[0].corners[i].x=mol->unitcell[mol->nhist].corners[i].x;
          mol->unitcell[0].corners[i].y=mol->unitcell[mol->nhist].corners[i].y;
          mol->unitcell[0].corners[i].z=mol->unitcell[mol->nhist].corners[i].z;
        }
      mol->unitcell[0].a=mol->unitcell[mol->nhist].a;
      mol->unitcell[0].b=mol->unitcell[mol->nhist].b;
      mol->unitcell[0].c=mol->unitcell[mol->nhist].c;
      mol->unitcell[0].alpha=mol->unitcell[mol->nhist].alpha;
      mol->unitcell[0].beta=mol->unitcell[mol->nhist].beta;
      mol->unitcell[0].gamma=mol->unitcell[mol->nhist].gamma;
        shiftUnitCell(mol);
        if (mol->unitcell) mol->thermoSettings&=~(TRANSLATION | ROTATION);
      }
      if (mol->history != NULL)
      {
        radfac1=1.0/radfac;
        for (i=0; i<mol->nhist*mol->na; i++)
        {
          mol->history[i].x*=radfac;
          mol->history[i].y*=radfac;
          mol->history[i].z*=radfac;
          mol->history[i].gx*=radfac1;
          mol->history[i].gy*=radfac1;
          mol->history[i].gz*=radfac1;
        }
      }
      if (mol->cnm != NULL)
      {
        for (i=0; i<mol->nmodes*mol->nmodes; i++)
          mol->cnm[i]*=radfac;
      }
      else
        mol->thermoSettings&=~VIBRATION;
/*    if (na != namax && namax != 0) setGeometry(FALSE); */
      imol=nmolecule;
      nmolecule++;
      return(TRUE);
    }
    else
      readNext=TRUE;

    if (readNext)
      word=fgets(line, MAXLENLINE, file);
    else
      word=line;
  }   /* end of "while (word != ..." */
  if (mol->na == 0)
  {
    word=getStringResource(topShell, "noCoordinates");
    sprintf(label, word, filename);
    messgb(topShell, 3, label, buttons1, 1);
    return(FALSE);
  }
  return(TRUE);
}

double *readMOs(struct MOLECULE *mol, FILE *file, int spin, int *ncoeff,
                char *line)
{
  GLdouble energy;
  double *cmo=NULL, *oldcmo, *coeff, e, occ;
  size_t mno=100, mnc=2000;
  int n=0;
  char sym[4], *word, *w;
  long offset;
  register int i, j, k, l, m;

  if (mol->orbitals == NULL)
  {
    mol->orbitals=(struct ORBITAL *)getmem(mno, sizeof(struct ORBITAL));
    cmo=(double *)getmem(mnc, sizeof(double));
    *ncoeff=i=j=0;
  }
  else
  {
    mno=mol->nbasfu+100;
    mol->orbitals=(struct ORBITAL *)expmem((void *)mol->orbitals, mno,
                                           sizeof(struct ORBITAL));
    mnc=mol->nbasfu*mol->nbasfu+mnc,
    cmo=(double *)expmem((void *)cmo, mnc, sizeof(double));
    i=j=mol->nbasfu;
  }
  while (TRUE)
  {
    fgets(line, MAXLENLINE, file);
    if (line[0] == '$') break;
    if ((word=strtok(line, " \t")) != NULL)
    {
      word=strtok(NULL, " \t");
      strcpy(mol->orbitals[i].symmetry, word);
      word=strtok(NULL, "=");
      word=strtok(NULL, " \t");
      mol->orbitals[i].energy=atof(word);
      mol->orbitals[i].spin=spin;
      mol->orbitals[i].occupation=0.0;
      word=strtok(NULL, "=");
      word=strtok(NULL, " \t");
      n=atoi(word);
      (*ncoeff)+=n;
    }
    if (n != 0)
    {
      mol->orbitals[i].coeff=&cmo[j];
      do
      {
        fgets(line, MAXLENLINE, file);
        if ((word=strrchr(line, '\n')) != NULL) *word='\0';
        if ((word=strtok(line, " \t")) != NULL)
        {
          cmo[j]=atof(word);
          n--;
          if (++j >= mnc)
          {
            mnc+=2000;
            oldcmo=cmo;
            cmo=(double *)expmem((void *)cmo, mnc, sizeof(double));
            if (oldcmo != cmo)
            {
              offset=(char *)cmo-(char *)oldcmo;
              for (l=0; l<=i; l++)
              {
                w=(char *)mol->orbitals[l].coeff+offset;
                mol->orbitals[l].coeff=(double *)w;
              }
            }
          }
        }
        while ((word=strtok(NULL, " \t")) != NULL)
        {
          cmo[j]=atof(word);
          n--;
          if (++j >= mnc)
          {
            mnc+=2000;
            oldcmo=cmo;
            cmo=(double *)expmem((void *)cmo, mnc, sizeof(double));
            if (oldcmo != cmo)
            {
              offset=(char *)cmo-(char *)oldcmo;
              for (l=0; l<=i; l++)
              {
                w=(char *)mol->orbitals[l].coeff+offset;
                mol->orbitals[l].coeff=(double *)w;
              }
            }
          }
        }
      } while (n > 0);
    }
    if (++i >= mno)
    {
      mno+=100;
      mol->orbitals=(struct ORBITAL *)expmem((void *)mol->orbitals, mno,
                                             sizeof(struct ORBITAL));
    }
  }  /* end of "while (TRUE)" */
  mol->nbasfu=i;
  if (mol->nbasfu > 0)
    mol->orbitals=(struct ORBITAL *)expmem((void *)mol->orbitals, (size_t)mol->nbasfu,
                                           sizeof(struct ORBITAL));
  oldcmo=cmo;
  if (j > 0)
  {
    cmo=(double *)expmem((void *)cmo, (size_t)j, sizeof(double));
    if (oldcmo != cmo)
    {
      k=(char *)cmo-(char *)oldcmo;
      for (l=0; l<i; l++)
      {
        m=(int)mol->orbitals[l].coeff+k;
        mol->orbitals[l].coeff=(double *)m;
      }
    }
  }
  for (i=0; i<mol->nbasfu; i++)
  {
    e=mol->orbitals[i].energy;
    k=i;
    for (j=i+1; j<mol->nbasfu; j++)
    {
      if (mol->orbitals[j].energy < e)
      {
        e=mol->orbitals[j].energy;
        k=j;
      }
    }
    if (k != i)
    {
      energy=mol->orbitals[k].energy;
      occ=mol->orbitals[k].occupation;
      coeff=mol->orbitals[k].coeff;
      l=mol->orbitals[k].spin;
      strcpy(sym, mol->orbitals[k].symmetry);
      mol->orbitals[k].energy=mol->orbitals[i].energy;
      mol->orbitals[k].occupation=mol->orbitals[i].occupation;
      mol->orbitals[k].coeff=mol->orbitals[i].coeff;
      mol->orbitals[k].spin=mol->orbitals[i].spin;
      strcpy(mol->orbitals[k].symmetry, mol->orbitals[i].symmetry);
      mol->orbitals[i].energy=energy;
      mol->orbitals[i].occupation=occ;
      mol->orbitals[i].coeff=coeff;
      mol->orbitals[i].spin=l;
      strcpy(mol->orbitals[i].symmetry, sym);
    }
  }
  return(cmo);
}

int readOccupation(struct MOLECULE *mol, FILE *file, int spin, char *line)
{
  int i, j, *help;
  char *word, *w;

  if (mol->nbasfu == 0)
  {
    word=fgets(line, MAXLENLINE, file);
    return(FALSE);
  }
  help=(int *)getmem(mol->nbasfu, sizeof(int));
  while (fgets(line, MAXLENLINE, file) != NULL)
  {
    if (strchr(line, '$')) break;
    word=strtok(line, " \t");
    w=strtok(NULL, " \t");
    i=getNumber(w, help);
    for (j=0; j<mol->nbasfu; j++)
    {
      if (!strcmp(mol->orbitals[j].symmetry, word) && mol->orbitals[j].spin == spin)
      {
        if (spin == ALPHAANDBETA)
          mol->orbitals[j].occupation=2.0;
        else
          mol->orbitals[j].occupation=1.0;
        if (--i == 0) break;
      }
    }
  }
  fremem((void **)&help);
  return(TRUE);
}

void relint(struct MOLECULE *mol)
{
  double wimax=0.0, rmimax=0.0;
  register int i;

/* Compute relative intensities (0 - 100 %) */
  for (i=0; i<mol->nmodes; i++)
  {
    if (mol->normal_modes[i].ir_intensity < 0.0)
      mol->normal_modes[i].ir_intensity=0.0;
    if (mol->normal_modes[i].raman_intensity < 0.0)
      mol->normal_modes[i].raman_intensity=0.0;
    wimax=MAX(wimax, mol->normal_modes[i].ir_intensity);
    rmimax=MAX(rmimax, mol->normal_modes[i].raman_intensity);
    if (mol->normal_modes[i].wavenumber < 0.0) mol->imag++;
  }
  if (wimax > 0.0)
    wimax=100./wimax;
  else
    wimax=1.0;
  if (rmimax > 0.0)
    rmimax=100./rmimax;
  else
    rmimax=1.0;
  for (i=0; i<mol->nmodes; i++)
  {
    mol->normal_modes[i].rel_ir_intensity=mol->normal_modes[i].ir_intensity*wimax;
    mol->normal_modes[i].rel_raman_intensity=mol->normal_modes[i].raman_intensity*rmimax;
    mol->normal_modes[i].rel_ins_intensity=0.0;
  }
}

int getNumber(char *str, int *array)
{
  register int i=0, j, cont=FALSE;
  register char *start, *p;

  start=p=str;
  while (*p)
  {
    if (*p == '-')
    {
      cont=TRUE;
      *p='\0';
      array[i++]=atoi(start);
      start=p+1;
    }
    else if (*p == ',')
    {
      *p='\0';
      if (cont)
      {
        for (j=array[i-1]+1; j<=atoi(start); j++)
          array[i++]=j;
        cont=FALSE;
      }
      else
        array[i++]=atoi(start);
      start=p+1;
    }
    p++;
  }
  if (cont)
  {
    for (j=array[i-1]+1; j<=atoi(start); j++)
      array[i++]=j;
  }
  else
    array[i++]=atoi(start);
  return(i);
}

int assignBasisset(struct MOLECULE *mol, int shell, char *label)
{
  int found=FALSE;
  register char *p;
  register int i;

  p=label+strlen(label)-1;
  while (*p == '\n' || *p == ' ') p--;
  *(p+1)='\0';
  p=label;
  while (*p == ' ') p++;
  for (i=0; i<mol->na; i++)
  {
    if (!strcmp(mol->atoms[i].basisname, p))
    {
      mol->atoms[i].basis=&mol->basisset[shell];
      *(mol->atoms[i].basisname)='\0';
      found=TRUE;
    }
  }
  return(found);
}

int moveBasissetData(struct MOLECULE *mol, struct BASISSET *old, int shell, int notify)
{
  static struct PushButtonRow buttons[] = {{"continue", GetMessageBoxButton, (XtPointer)0, NULL}};
  size_t k;
  char *missing=NULL, number[10], *str;
  int ierror, count=0;
  register long offset;
  register int i, j;
  register char *p;

  ierror=FALSE;
  offset=(char *)(mol->basisset)-(char *)old;
  k=0;
  for (j=0; j<mol->na; j++)
  {
    if (mol->atoms[j].basis != NULL)
    {
      p=(char *)mol->atoms[j].basis+offset;
      mol->atoms[j].basis=(struct BASISSET *)p;
    }
    else
    {
      if (notify)
      {
        k+=12;
        if (j < 1000000)
          sprintf(number, "%s(%d)", mol->atoms[j].name, j+1);
        else
          sprintf(number, "%s", mol->atoms[j].name);
        if (missing == NULL)
        {
          missing=(char *)getmem(k, sizeof(char));
          count=0;
        }
        else
        {
          missing=(char *)expmem((void *)missing, k, sizeof(char));
          strcat(missing, ", ");
          count+=12;
          if (count > 60)
          {
            count=0;
            missing[strlen(missing)-1]='\n';
          }
        }
        strcat(missing, number);
      }
    }
  }
  if (k > 0 && showWarning)
  {
    ierror=TRUE;
    p=getStringResource(topShell, "missingBasisset");
    if (!p) p="The basis set for the atom(s) %s has not been specified.";
    str=(char *)getmem((size_t)(strlen(p)+strlen(missing)+1), sizeof(char));
    sprintf(str, p, missing);
    messgb(topShell, 1, str, buttons, 1);
    fremem((void **)&str);
    fremem((void **)&missing);
  }
  for (i=0; i<shell; i++)
  {
    mol->basisset[i].first=old[i].first;
    mol->basisset[i].nprim=old[i].nprim;
    mol->basisset[i].ang=old[i].ang;
    mol->basisset[i].next=old[i].next;
    if (mol->basisset[i].next != NULL)
    {
      p=(char *)(mol->basisset[i].next)+offset;
      mol->basisset[i].next=(struct BASISSET *)p;
    }
  }
  return(ierror);
}

void changeExponentPointers(struct BASISSET *basisset, long offset, int shell)
{
  register int i;

  for (i=0; i<=shell; i++)
    basisset[i].first+=offset;
}

void sortBasis(struct MOLECULE *mol)
{
  struct BASISSET *shell1, *shell2, *sorted;
  int type;
  register int i, j, k;

  for (i=0; i<mol->na; i++)
  {
    j=0;
    shell1=mol->atoms[i].basis;
    while (shell1 != NULL)
    {
      j++;
      shell1=shell1->next;
    }
    sorted=(struct BASISSET *)getmem(j, sizeof(struct BASISSET));
    k=0;
    shell1=mol->atoms[i].basis;
    while (shell1 != NULL)
    {
      type=shell1->ang;
      if (type != -1)
      {
        shell2=shell1;
        while (shell2 != NULL)
        {
          if (shell2->ang == type)
          {
            sorted[k].first=shell2->first;
            sorted[k].ang=shell2->ang;
            sorted[k++].nprim=shell2->nprim;
            shell2->ang=(-1);
          }
          shell2=shell2->next;
        }
      }
      shell1=shell1->next;
    }
    k=0;
    shell1=mol->atoms[i].basis;
    while (shell1 != NULL)
    {
      shell1->first=sorted[k].first;
      shell1->ang=sorted[k].ang;
      shell1->nprim=sorted[k++].nprim;
      shell1=shell1->next;
    }
    free(sorted);
  }
}

int compar(const void *first, const void *second)
{
  struct BASISSET *f=(struct BASISSET *)first, *s=(struct BASISSET *)second;

  return(f->ang - s->ang);
}

struct SAVE *saveGeometry(struct MOLECULE *mol)
{
  struct SAVE *coord;
  register int i;

  coord=(struct SAVE*)getmem(mol->na, sizeof(struct SAVE));
  for (i=0; i<mol->na; i++)
  {
    coord[i].x=mol->atoms[i].x;
    coord[i].y=mol->atoms[i].y;
    coord[i].z=mol->atoms[i].z;
    strcpy(coord[i].name, mol->atoms[i].name);
  }
  return(coord);
}

void wrbas(struct MOLECULE *mol)
{                                    /* dump basis set for debugging */
  struct BASISSET *shell;
  int *done;
  register int i, j;
  register char ang;
  register double *p;

  done=(int *)getmem((size_t)mol->na, sizeof(int));

  for (i=0; i<mol->na; i++)
  {
    if (done[i]) continue;
    printf("Atom(s) %d", i+1);
    done[i]=TRUE;
    shell=mol->atoms[i].basis;
    for (j=i+1; j<mol->na; j++)
    {
      if (mol->atoms[j].basis == shell)
      {
        printf(", %d", j+1);
        done[j]=TRUE;
      }
    }
    printf("\n");
    do
    {
      switch (shell->ang)
      {
        case 0:  ang='s';
                 break;
        case 1:  ang='p';
                 break;
        case 2:  ang='d';
                 break;
        default: ang=shell->ang+'c';
                 break;
      }
      printf("%4d  %c\n", shell->nprim, ang);
      p=shell->first;
      for (j=0; j<shell->nprim; j++)
      {
        printf("%f  %f\n", *p, *(p+1));
        p+=2;
      }
      shell=shell->next;
    } while (shell != NULL);
  }
  fremem((void **)&done);
}

struct MOLECULE *initMolecule(int mna)
{
  PyMoleculeSpecObject *newMolecule;
  struct MOLECULE *mol;
  register int i, j;

  if (molecules)
    molecules=(struct MOLECULE *)expmem((void *)molecules, (size_t)(nmolecule+1),
                                        sizeof(struct MOLECULE));
  else
    molecules=(struct MOLECULE *)getmem((size_t)1, sizeof(struct MOLECULE));
  mol=&molecules[nmolecule];
  if ((newMolecule=molecule_new()) != NULL)
  {
    newMolecule->moleculeID=nmolecule;
    mol->pyObject=newMolecule;
  }
  mol->atoms=(struct ATOM *)getmem(mna, sizeof(struct ATOM));
  mol->bonds=NULL;
  mol->unitcell=NULL;
  mol->normal_modes=NULL;
  mol->history=NULL;
  mol->internals=NULL;
  mol->optimization=NULL;
  mol->orbitals=NULL;
  mol->gridObjects=NULL;
  mol->basisset=NULL;
  mol->cnm=NULL;
  mol->exponents=NULL;
  mol->coord=NULL;
  setBond(mol);
  mol->addedBonds=NULL;
  mol->deletedBonds=NULL;
  mol->bondShift=0.0;
  mol->transx=0.0;
  mol->transy=0.0;
  mol->transz=0.0;
  mol->rotConstants[0]=0.0;
  mol->rotConstants[1]=0.0;
  mol->rotConstants[2]=0.0;
  for (i=0; i<3;i++)
  {
    for (j=0; j<3; j++)
      mol->tinert[i][j]=0.0;
  }
/*mol->cellFactor[0]=1.0;
  mol->cellFactor[1]=1.0;
  mol->cellFactor[2]=1.0; */
  mol->na=0;
  mol->nb=0;
  mol->nbAdded=0;
  mol->nbDeleted=0;
  mol->nmodes=0;
  mol->nhist=0;
  mol->ninternal=0;
  mol->nbasfu=0;
  mol->ngridobjects=0;
  mol->pured=TRUE;
  mol->puref=TRUE;
  mol->pureg=TRUE;
  mol->reaction=NOTHING;
  mol->stoichioNumber=0;
  mol->imag=0;
  mol->mode=(-1);
  mol->imo=(-1);
  mol->imosave=(-1);
  mol->ibasfu=(-1);
  mol->needMoloch=FALSE;
  mol->gaussian=FALSE;
  mol->thermoSettings=TRANSLATION | PV | ROTATION | VIBRATION;
/*mol->showMiller=FALSE;
  mol->miller[0]=0; mol->miller[1]=mol->miller[2]=1;
  mol->nmiller=0; */
  mol->title[0]='\0';
/*sprintf(mol->title, "Molecule %d", nmolecule+1); */
  sprintf(mol->pgroup, "c1");
  transObject=(float *)expmem((void *)transObject, nmolecule+MOLECULES+1, 3*sizeof(float));
  i=3*(nmolecule+MOLECULES);
  transObject[i]=transObject[i+1]=transObject[i+2]=0.0;
  rotObject=(float *)expmem((void *)rotObject, nmolecule+MOLECULES+1, 4*sizeof(float));
  i=4*(nmolecule+MOLECULES);
  rotObject[i]=rotObject[i+1]=rotObject[i+2]=0.0;
  rotObject[i+3]=1.0;
  return(mol);
}

void deleteMolecule(Widget w, caddr_t interactive, XmAnyCallbackStruct *callbackData)
{
  static struct PushButtonRow buttons[] = {{"ok", GetMessageBoxButton, (XtPointer)TRUE, NULL},
                                           {"cancel", GetMessageBoxButton, (XtPointer)FALSE, NULL}};
  int first, last;
  char *word, line[MAXLENLINE];
  register struct MOLECULE *mol;
  register int i, j;

  first=windows[VIEWER].set;
  last=first+1;
  if ((int)interactive)
  {
    if (first == -1)
    {
      first=0;
      last=nmolecule;
      word=getStringResource(topShell, "deleteAll");
      if (!messgb(topShell, 2, word, buttons, 2)) return;
    }
    else
    {
      word=getStringResource(topShell, "deleteOne");
      sprintf(line, word, molecules[first].title);
      if (!messgb(topShell, 2, line, buttons, 2)) return;
    }
  }

  for (i=first; i<last; i++)
  {
    mol=&molecules[first];
    Py_DECREF(mol->pyObject);
    fremem((void **)&mol->atoms);
    if (mol->bonds != NULL)        fremem((void **)&(mol->bonds));
    if (mol->unitcell != NULL)     fremem((void **)&(mol->unitcell));
    if (mol->normal_modes != NULL) fremem((void **)&(mol->normal_modes));
    if (mol->history != NULL)      fremem((void **)&(mol->history));
    if (mol->internals != NULL)    fremem((void **)&(mol->internals));
    if (mol->optimization != NULL) fremem((void **)&(mol->optimization));
    if (mol->orbitals != NULL)     fremem((void **)&(mol->orbitals));
    if (mol->gridObjects != NULL)  fremem((void **)&(mol->gridObjects));
    if (mol->basisset != NULL)     fremem((void **)&(mol->basisset));
    if (mol->coord != NULL)        fremem((void **)&(mol->coord));
    if (mol->cnm != NULL)          fremem((void **)&(mol->cnm));
    if (mol->exponents != NULL)    fremem((void **)&(mol->exponents));
    if (mol->addedBonds != NULL)   fremem((void **)&(mol->addedBonds));
    if (mol->deletedBonds != NULL) fremem((void **)&(mol->deletedBonds));
/*  PyMem_DEL(mol->pyObject);*/

    if (first != nmolecule-1)
    {
      for (j=first+1; j<nmolecule; j++)
      {
        molecules[j-1]=molecules[j];
        (molecules[j-1].pyObject)->moleculeID=j-1;
        transObject[3*(j+MOLECULES-1)]=transObject[3*(j+MOLECULES)];
        transObject[3*(j+MOLECULES-1)+1]=transObject[3*(j+MOLECULES)+1];
        transObject[3*(j+MOLECULES-1)+2]=transObject[3*(j+MOLECULES)+2];
        rotObject[4*(j+MOLECULES-1)]=rotObject[4*(j+MOLECULES)];
        rotObject[4*(j+MOLECULES-1)+1]=rotObject[4*(j+MOLECULES)+1];
        rotObject[4*(j+MOLECULES-1)+2]=rotObject[4*(j+MOLECULES)+2];
        rotObject[4*(j+MOLECULES-1)+3]=rotObject[4*(j+MOLECULES)+3];
      }
    }
    nmolecule--;
    if (nmolecule)
      molecules=(struct MOLECULE *)expmem((void *)molecules, (size_t)nmolecule,
                                          sizeof(struct MOLECULE));
    else
      fremem((void **)&molecules);
    transObject=(float *)expmem((void *)transObject, (size_t)(nmolecule+MOLECULES),
                                3*sizeof(float));
    rotObject=(float *)expmem((void *)rotObject, (size_t)(nmolecule+MOLECULES),
                              4*sizeof(float));
    windows[VIEWER].set=nmolecule-1;
  }
  windows[VIEWER].menu=makeViewerMenu(windows[VIEWER].widget);

  if (nmolecule < 2) setMenuItem(VIEWER_SELECT, False);
  if (nmolecule == 0)
  {
    setMenu((struct MOLECULE *)NULL);
    setMenuItem(VIEWER1_SAVE, False);
    setMenuItem(VIEWER1_DELETE, False);
    setMenuItem(VIEWER1_BUILD, False);
    for (i=VIEWER_GEOMETRY; i<VIEWER_SEP3; i++)
      setMenuItem(i, False);
    setMenuItem(VIEWER_HARDCOPY, False);
    word=getStringResource(topShell, "title");
    setWindowTitle(windows[VIEWER].widget, word);
    XtRemoveAllCallbacks(windows[VIEWER].widget, XmNexposeCallback);
    XtAddCallback(windows[VIEWER].widget, XmNexposeCallback, (XtCallbackProc)showTitle, NULL);
    if (windows[SPECTRUM].widget) quitSpectrum((Widget)0, (caddr_t)0,
                                               (XmAnyCallbackStruct *)0);
    if (windows[HISTORY].widget) quitHistory((Widget)0, (caddr_t)0,
                                             (XmAnyCallbackStruct *)0);
    if (windows[MO].widget) quitMODiagram((Widget)0, (caddr_t)0,
                                          (XmAnyCallbackStruct *)0);
  }
  else
    setWindowTitle(windows[VIEWER].widget, molecules[windows[VIEWER].set].title);
  redraw(VIEWER);
}

void reformat(char *source, char *destination)
{
  static struct PushButtonRow buttons[] = {{"exit", GetMessageBoxButton,
                                            (XtPointer)0, NULL}};
  FILE *s, *d;
  char line[MAXLENLINE], *word;

  s=fopen(source, "r");
  if ((d=fopen(destination, "w")) == NULL)
  {
    word=getStringResource(topShell, "cannotOpen");
    sprintf(line, word, destination);
    messgb(topShell, 3, line, buttons, 1);
    exit(-1);
  }
  while (fgets(line, MAXLENLINE, s))
  {
    if ((word=strrchr(line, '\r')) != NULL)
      *word='\0';
    fprintf(d, "%s\n", line);
  }
  fclose(d);
  fclose(s);
}