/**************************************************************************/
/**************************************************************************/
/**************************   "steric"   **********************************/
/**************************************************************************/
/*************     Program to calculate ligand cone    ********************/
/*************     angles as a measure of steric size  ********************/
/**************************************************************************/
/**************************************************************************/

/**************************************************************************/
/******************        Main Function         **************************/
/**************************************************************************/
/******************        This module is        **************************/
/******************      system independant      **************************/
/**************************************************************************/

#include 
#include 
#include 

#include "sterdefn.h" /* declaration of structures, functions and globals */
#include "stercomm.h" /* definitions of all menu command options          */
#include "crystal.h"       /* functions for all file manipulations        */
#include "stercalc.h" /* main calculation routines                        */
#include "stermem.h"  /* dynamic memory management module                 */
#include "sterfile.h" /* file input and output functions                  */
#include "stertext.h" /* all text functions for text mode                 */
#include "steraid.h"  /* additional functions needed                      */
#include "stergrap.h" /* graph plotting funtions                          */
#include "stergrp.h"  /* group manipulation functions                     */

FILE *output;         /* pointer to text output                           */

/**************************************************************************/
/**********************   file operations  ********************************/
/**************************************************************************/

Mol *Load_Molecules(char *line, Mol *molecule, Set *set)
{
  char *args[MAXARG];
  char file[LINELEN];
  int num=0,count=0;
  Mol *new=NULL;

  if((num=Get_Arguements(line,args))>0)
  for(count=0;countFname);
  }
}

/**************************************************************************/

void Open_Shell(void)
{
  Out_Message("Entering a unix csh, if available",O_NEWLN);
  Out_Message("All normal csh commands should be available",O_NEWLN);
  Out_Message("Type \"exit\" to return to steric",O_NEWLN);
  system("csh");
}

/**************************************************************************/
/**********************  Molecule operations  *****************************/
/**************************************************************************/

Mol *Test_Close_Current_Molecule(Mol *M, char *argline)
{
  int a,n;
  char *args[MAXARG];
  int num=0;
  Grps *group=NULL;
  char line[LINELEN];
  char inline[LINELEN];
  char outline[LINELEN];

  if(M==NULL) return(NULL);
  if(strlen(argline)<1) return(Close_Current_Molecule(M));
  if(M->groups==NULL)
  {
    Error_Message(E_NOGRPS,"Test Close Current Molecule");
    return(M);
  }

/* load specified group names into memory */

  strcpy(inline,argline);
  if((num=Get_Group_Names(inline,line,args))<1) return(M);

/* loop through all molecule groups, comparing to loaded groups           */

  for(a=0,n=0;ngroups);group!=NULL;group=group->next)
    {
      if(stricmp(group->name,args[n])==0)
      {
        sprintf(outline,"Molecule %s has group named %s",M->name,args[n]);
        Out_Message(outline,O_NEWLN);
        a++;
        break;
      }
    }
    if(group==NULL)
    {
      sprintf(outline,"Molecule %s doesn't have group %s",M->name,args[n]);
      Out_Message(outline,O_NEWLN);
    }
  }
  if(a!=num)
  {
    sprintf(outline,"Molecule %s doesn't have all required groups",M->name);
    Out_Message(outline,O_NEWLN);
    M=Close_Current_Molecule(M);
  }
  return(M);
}

/**************************************************************************/

Mol *Goto_Particular_Molecule(Mol *molecule, char *line, Set *set)
{
  char *args[MAXARG];
  char string[LINELEN];
  int num=0;
  Mol *M=NULL;
  int mi=0;

  if(molecule==NULL) return(NULL);
  if(!strlen(line))
  {
    sprintf(line,"Enter number or label of molecule to make current:");
    Get_Input_Line(line,set->input);
  }
  if((strlen(line))&&((num=Get_Arguements(line,args))>0))
  {
    strcpy(string,args[0]);
    if(((M=Goto_Molecule_Name(molecule,string))==NULL)
     &&((mi=Int_String(string))==0)) Error_Message(E_UNKNOW,"Kill Atom");
    else
    {
      if((M==NULL)&&(mi)) M=Goto_Molecule(molecule,mi);
    }
    if(M!=NULL) molecule=M;
  }
  return(molecule);
}

/**************************************************************************/

void Kill_Atoms(Mol *M, char *line, FILE *input)
{
  char *args[MAXARG];
  char string[LINELEN];
  int num=0,count=0;
  Atms *at=NULL;
  int ati=0;

  if(!strlen(line))
  {
    sprintf(line,"Enter numbers or labels of atoms to kill:");
    Get_Input_Line(line,input);
  }
  if(strlen(line))
  {
    if((num=Get_Arguements(line,args))>0)
    for(count=0;countatoms,string))==NULL)
       &&((ati=Int_String(string))==0)) Error_Message(E_UNKNOW,"Kill Atom");
      else
      {
        if(at) at->stat=at->stat&(!MAIN_BIT);
        else if(ati)
        {
          at=Goto_Atom(M->atoms,ati);
          at->stat=at->stat&(FULL_BIT^MAIN_BIT);
        }
      }
    }
  }
}

/**************************************************************************/

void Kill_Types(Mol *M, char *line, FILE *input)
{
  char *args[MAXARG];
  char string[LINELEN];
  int num=0,count=0;
  Atms *at=NULL;

  if(!strlen(line))
  {
    sprintf(line,"Enter labels of atom types to kill:");
    Get_Input_Line(line,input);
  }
  if(strlen(line))
  {
    if((num=Get_Arguements(line,args))>0)
    for(count=0;countatoms);at!=NULL;at=at->next)
      {
        if(strcmp(string,at->type)==0) at->stat=at->stat&(FULL_BIT^MAIN_BIT);
      }
    }
  }
}

/**************************************************************************/

void Kill_Groups_Name(Mol *M, char *line, FILE *input)
{
  char *args[MAXARG];
  char string[LINELEN];
  int num=0,count=0;
  Atms *at=NULL;
  Grps *gr=NULL;
  int gri=0;

  if(!strlen(line))
  {
    sprintf(line,"Enter numbers or labels of groups to kill:");
    Get_Input_Line(line,input);
  }
  if(strlen(line))
  {
    if((num=Get_Arguements(line,args))>0)
    for(count=0;countgroups,string))==NULL)) Error_Message(E_UNKNOW,"Kill Groups Name");
      else
      {
        if(gr==NULL) gr=Goto_Group(M->groups,gri);
        if(gr!=NULL)
        {
          gri=Get_Group_Number(gr);
          for(at=First_Atom(M->atoms);at!=NULL;at=at->next)
          {
            if((at!=NULL)&&(at->group==gri)) at->stat=at->stat&(!MAIN_BIT);
          }
        }
      }
    }
  }
}

/**************************************************************************/

int Kill_Groups_Type(Mol *M, char *line, FILE *input)
{
  char *args[MAXARG];
  char string[LINELEN];
  int num=0,count=0;

  if(M->groups==NULL)
  {
    Error_Message(E_NOGRPS,"Kill Groups Type");
    return(0);
  }
  if(!strlen(line))
  {
    sprintf(line,"Enter labels of group types to kill:");
    Get_Input_Line(line,input);
  }
  if(strlen(line))
  {
    if((num=Get_Arguements(line,args))>0)
    for(count=0;countorigin!=NULL)
  {
    sprintf(line,"Only origin defining atom will remain excluded");
    Out_Message(line,O_NEWLN);
  }
  for(at=First_Atom(M->atoms);at!=NULL;at=at->next)
  {
    if(at!=M->origin) at->stat=at->stat|MAIN_BIT;
  }
}

/**************************************************************************/

void Kill_Molecule_Data(char *args, char *line, Mol *M, Set *set)
{
  switch(args[0])
  {
    case ATMS : Kill_Atoms(M,line,set->input); break;
    case TYPS : Kill_Types(M,line,set->input); break;
    case GRPS : Kill_Groups_Type(M,line,set->input); break;
    case GRPN : Kill_Groups_Name(M,line,set->input); break;
    case REST : UnKill_All(M,line); break;
    default   : break;
  }
  if(M->origin!=NULL) Calculate_Parameters(M,M->origin->v,set);
}

/**************************************************************************/
/**********************    calculations    ********************************/
/**************************************************************************/

int Calculate_Steric_Value(char *args, Mol *molecule, Set *set)
{
  Ster *S=NULL;
  int i,len=strlen(args);
  for(i=0;ister,set);
    if(S!=NULL)
    {
      molecule->ster=S;
      molecule->multi=set->multi;
      Calc_Total_Steric(molecule,S,set,args[i]);
    }
  }
  return(i);
}

/**************************************************************************/

int Calculate_Steric_Conformer_Average(char *args, Mol *molecule, Set *set)
{
  Ster *S=NULL;
  int i,len=strlen(args);
  for(i=0;ister,set);
    if(S!=NULL)
    {
      molecule->ster=S;
      molecule->multi=set->multi;
      Calc_Conformer_Average(molecule,S,set,args[i]);
    }
  }
  return(i);
}

/**************************************************************************/

int Calculate_Steric_Profile(char *args, Mol *molecule, Set *set)
{
  Ster *S=NULL;
  int i,len=strlen(args);
  for(i=0;ister,set);
    if(S!=NULL)
    {
      molecule->ster=S;
      molecule->multi=set->multi;
      if(S->type==S_ANGU)
        Calc_Angular_Profile(molecule,S,set);
      else
        Calc_Radial_Profile(molecule,S,set,args[i]);
    }
  }
  return(i);
}

/**************************************************************************/

int Calculate_Volumes(char *args, Mol *molecule, Set *set, char *line)
{
  switch(args[0])
  {
    case FREV : Calculate_Free_Unit_Cell_Volume(molecule,line,set); break;
    case GRPV : Calculate_Group_Volume(molecule,line,set,GRPV); break;
    case CAVV : Calculate_Group_Volume(molecule,line,set,CAVV); break;
    default   : break;
  }
  return(1);
}

/**************************************************************************/

int Calculate_Areas(char *args, Mol *molecule, Set *set)
{
  switch(args[0])
  {
    case THET : Calculate_Areas_OneD(molecule,set,args[0]); break;
    case PHIP : Calculate_Areas_OneD(molecule,set,args[0]); break;
    case MAPP : Calculate_Areas_TwoD(molecule,set); break;
    default   : args[0]=PROJ; args[1]=0; Calculate_Steric_Value(args,molecule,set); break;
  }
  return(1);
}

/**************************************************************************/

int Calculate_Group_Shell(char *args, Mol *M, Set *set, char *argline)
{
  Ster *S=NULL;
  int i,len=strlen(args);
  int a,n,count,gnum;
  Atms *A=NULL;
  char *gargs[MAXARG];
  int num=0;
  Grps *group=NULL;
  char line[LINELEN];
  char inline[LINELEN];
  double smooth=0.0;
  unsigned Mmode;

  if(M==NULL) return(0);
  if(strlen(argline)<1) return(0);
  if(M->groups==NULL)
  {
    Error_Message(E_NOGRPS,"Calculate Group Shell");
    return(0);
  }

  strcpy(inline,argline);
  if((num=Get_Group_Names(inline,line,gargs))<1) return(0);

  Mmode=set->mode;
  set->mode&=FULL_BIT^RSET_BIT;
  set->mode&=FULL_BIT^RDEF_BIT;
  for(i=0;ister,set);
    if(S->type==S_ANGU)
    {
      M->ster=Close_Steric(S);
      break;
    }
    if(S!=NULL)
    {
      M->ster=S;
      M->multi=set->multi;
      for(n=0;ngroups);group!=NULL;group=group->next)
        {
          if(stricmp(group->name,gargs[n])!=0) continue;
          if(group->main==NULL) continue;
          gnum=Get_Group_Number(group);
          M->origin=group->main;
          if(!AlmostZero(group->volrad))
          {
            M->minR=0.0;
            M->maxR=group->volrad*1.5;
            S->min=M->minR;
            S->max=M->maxR;
          }
          group->mode|=RADV_BIT;

          for(A=First_Atom(M->atoms);A!=NULL;A=A->next)
          {
            if(A->group==gnum) A->stat&=FULL_BIT^MAIN_BIT;
          }
          if(M->origin) Get_Positions(M,M->origin->v,set);
          Calc_Radial_Profile(M,S,set,args[i]);
          for(A=First_Atom(M->atoms);A!=NULL;A=A->next)
          {
            if((A!=M->origin)&&(A->group==gnum)) A->stat&=MAIN_BIT;
          }
          group->volrad=S->peak_R;
          for(a=set->smooth;asize-set->smooth;a++)
          {
            smooth=0.0;
            for(count=-set->smooth;count<=set->smooth;count++)
              smooth+=S->val[a+count];
            smooth/=(double)(1+2*set->smooth);
            if(smooth>=set->vrcut) 
            {
              group->volrad=(double)a*(S->max-S->min)/((double)S->size);
              break;
            }
          }
          group->radvol=4.0/3.0*PI*group->volrad*group->volrad*group->volrad;
          sprintf(inline,"%s group %d %s has shell at %f using %f cutoff, sphere volume = %f"
                        ,M->name,gnum,group->name,group->volrad,set->vrcut,group->radvol);
          Out_Message(inline,O_NEWLN);
        }
      }
    }
  }
  set->mode=Mmode;
  return(i);
}

/**************************************************************************/

int Calculate_Crystal_Volumes(char *args, Mol *molecule, Set *set, char *line)
{
  switch(args[0])
  {
    case FREV : Calculate_Free_Unit_Cell_Volume(molecule,line,set); break;
    case GRPV : Calculate_Group_Volume(molecule,line,set,GRPV); break;
    case CAVV : Calculate_Group_Volume(molecule,line,set,CAVV); break;
    default   : break;
  }
  return(1);
}

/**************************************************************************/

int Calculate_Crystal_Areas(char *args, Mol *molecule, Set *set, char *line)
{
  Ster *S=NULL;
  S=Which_Ster(PROJ,molecule->ster,set);
  if(S!=NULL)
  {
    molecule->ster=S;
    molecule->multi=set->multi;
    Calc_Total_Steric(molecule,S,set,PROJ);
  }
  return(1);
}

/**************************************************************************/
/********************** change data settings ******************************/
/**************************************************************************/

void Change_Profrange(Set *set)
{
  set->min=Change_Double("Enter minimum radial distance",set->min,set->input);
  set->max=Change_Double("Enter maximum radial distance",set->max,set->input);
  set->a_min=2.0*PI*Change_Double("Enter minimum profile angle (n*2pi)",set->a_min/(2.0*PI),set->input);
  set->a_max=2.0*PI*Change_Double("Enter maximum profile angle (n*2pi)",set->a_max/(2.0*PI),set->input);
}

/**************************************************************************/

void Change_Profarray_Size(Set *set)
{
  set->size=Change_Integer("Enter profile array size",set->size,set->input);
  set->t_size=Change_Integer("Enter area vs. theta profile array size",set->t_size,set->input);
  set->p_size=Change_Integer("Enter area vs. phi profile array size",set->p_size,set->input);
  set->n_size=Change_Integer("Enter number of steps in numerical solid angle calculation",set->n_size,set->input);
  set->c_size=Change_Integer("Enter number of angular steps in 3D profile plot",set->c_size,set->input);
  set->smooth=Change_Integer("Enter new profile smoothing window",set->smooth,set->input);
}

/**************************************************************************/

void Change_Bonding(Set *set)
{
  set->bond_minf=Change_Double("Enter fraction of covalent radius to be minimum bond length",set->bond_minf,set->input);
  set->bond_maxf=Change_Double("Enter fraction of covalent radius to be maximum bond length",set->bond_maxf,set->input);
}

/**************************************************************************/

void Change_Volume_Grid(Set *set)
{
  set->vx=Change_Integer("Enter volume x grid size",set->vx,set->input);
  set->vy=Change_Integer("Enter volume y grid size",set->vy,set->input);
  set->vz=Change_Integer("Enter volume z grid size",set->vz,set->input);
  set->tgap=Change_Integer("Enter MCV output time gap in seconds",set->tgap,set->input);
  set->maxmcv=Change_Integer("Enter MCV maximum number of grid points",set->maxmcv,set->input);
  set->minmcv=Change_Integer("Enter MCV minimum number of grid points",set->minmcv,set->input);
  set->avmmcv=Change_Integer("Enter MCV number of grid points at which to start averaging",set->minmcv,set->input);
}

/**************************************************************************/

void Change_Profile_Mode(Set *set)
{
  if(Get_Yes_No_Request("Allow general settings profile range to overide profile range?",set->mode&RSET_BIT,set->input))
    set->mode=set->mode|RSET_BIT;
  else if(Get_Yes_No_Request("Allow original molecule profile range to overide profile range?",set->mode&RDEF_BIT,set->input))
    set->mode=set->mode|(RDEF_BIT&(!RSET_BIT));
  else set->mode=(set->mode&(FULL_BIT^RSET_BIT))&(FULL_BIT^RDEF_BIT);
}

/**************************************************************************/

void Change_General_Mode(Set *set)
{
  Change_Profile_Mode(set);
  if(Get_Yes_No_Request("Use new G positioning in old leach solid angles?",set->mode&SCOR_BIT,set->input))
    set->mode=set->mode|MOVG_BIT;
  else set->mode=set->mode&(FULL_BIT^MOVG_BIT);
  if(Get_Yes_No_Request("Use faulty segment correction in old leach solid angles?",set->mode&SCOR_BIT,set->input))
    set->mode=set->mode|SCOR_BIT;
  else set->mode=set->mode&(FULL_BIT^SCOR_BIT);
  if(Get_Yes_No_Request("Use only non-bonded overlaps in vertex and solid angles of overlap?",set->mode&NBO_BIT,set->input))
    set->mode=set->mode|NBO_BIT;
  else set->mode=set->mode&(FULL_BIT^NBO_BIT);
  if(Get_Yes_No_Request("Use only overlaps between different groups in VAO and SAO?",set->mode&NGR_BIT,set->input))
    set->mode=set->mode|NGR_BIT;
  else set->mode=set->mode&(FULL_BIT^NGR_BIT);
  if(Get_Yes_No_Request("Output all cone contour data in numerical solid angle calculation?",set->mode&CONT_BIT,set->input))
    set->mode=set->mode|CONT_BIT;
  else set->mode=set->mode&(FULL_BIT^CONT_BIT);
  if(Get_Yes_No_Request("Use perspective mode in contour plots?",set->mode&PERS_BIT,set->input))
    set->mode=set->mode|PERS_BIT;
  else set->mode=set->mode&(FULL_BIT^PERS_BIT);
  if(Get_Yes_No_Request("Use van der Waals radii for all atoms?",set->mode&VDW_BIT,set->input))
    set->mode=set->mode|VDW_BIT;
  else set->mode=set->mode&(FULL_BIT^VDW_BIT);
  if(!(set->mode&VDW_BIT))
  {
    if(Get_Yes_No_Request("Use covalent radii for all atoms?",set->mode&CVL_BIT,set->input))
      set->mode=set->mode|CVL_BIT;
    else set->mode=set->mode&(FULL_BIT^CVL_BIT);
  }
  if(Get_Yes_No_Request("Use Monte Carlo approach to volume calculation?",set->mode&MTC_BIT,set->input))
    set->mode=set->mode|MTC_BIT;
  else set->mode=set->mode&(FULL_BIT^MTC_BIT);

  if(Get_Yes_No_Request("Allow dummy unit cell atoms to be visable?",set->modev&UCAV_BIT,set->input))
    set->modev=set->modev|UCAV_BIT;
  else set->modev=set->modev&(FULL_BIT^UCAV_BIT);
  if(Get_Yes_No_Request("Automatically box unit cell in cavity volume?",set->modev&ABOX_BIT,set->input))
    set->modev=set->modev|ABOX_BIT;
  else set->modev=set->modev&(FULL_BIT^ABOX_BIT);
  if(Get_Yes_No_Request("Automatically remove redundant atoms after box?",set->modev&REDU_BIT,set->input))
    set->modev=set->modev|REDU_BIT;
  else set->modev=set->modev&(FULL_BIT^REDU_BIT);
  if(Get_Yes_No_Request("Use atom name in duplication test?",set->modev&SNAM_BIT,set->input))
    set->modev=set->modev|SNAM_BIT;
  else set->modev=set->modev&(FULL_BIT^SNAM_BIT);
  if(Get_Yes_No_Request("View 2D coordinates in atom listing?",set->modev&TWOD_BIT,set->input))
    set->modev=set->modev|TWOD_BIT;
  else set->modev=set->modev&(FULL_BIT^TWOD_BIT);
  if(Get_Yes_No_Request("View fractional coordinates in atom listing?",set->modev&FRAX_BIT,set->input))
    set->modev=set->modev|FRAX_BIT;
  else set->modev=set->modev&(FULL_BIT^FRAX_BIT);
  if(Get_Yes_No_Request("View polar coordinates in atom listing?",set->modev&POLA_BIT,set->input))
    set->modev=set->modev|POLA_BIT;
  else set->modev=set->modev&(FULL_BIT^POLA_BIT);
  if(Get_Yes_No_Request("Fade usused atom in orthogonal plot?",set->modev&FADE_BIT,set->input))
    set->modev=set->modev|FADE_BIT;
  else set->modev=set->modev&(FULL_BIT^FADE_BIT);
}

/**************************************************************************/

void Change_Accuracy(Set *set)
{
  set->eps=Change_Double("Enter minimum error for Simpsons rule integration",set->eps,set->input);
  set->veps=Change_Double("Enter minimum error for Monte Carlo volume",set->veps,set->input);
  set->vrcut=4.0*PI*Change_Double("Enter new volrad profile cutoff (n*4pi)",set->vrcut/(4.0*PI),set->input);
}

/**************************************************************************/

void Change_Contour_Rotation(Set *set)
{
  set->rotation=2*PI*Change_Double("Enter new contour plot rotation angle about z (n*2pi)",set->rotation/(2*PI),set->input);
}

/**************************************************************************/

void Change_Order_Multi(Set *set)
{
  set->multi=Change_Integer("Enter new maximum overlap order to consider in Craig solid angle",set->multi,set->input);
}

/**************************************************************************/

void Change_Plot_Parameters(Set *set)
{
  set->x_ran=Change_Integer("Enter new integer x range",set->x_ran,set->input);
  set->y_ran=Change_Integer("Enter new integer y range",set->y_ran,set->input);
}

/**************************************************************************/

void Change_General_Settings(char *args, Set *set)
{
  switch(args[0])
  {
    case ORDM : Change_Order_Multi(set);
    case RANG : Change_Profrange(set); break;
    case SIZE : Change_Profarray_Size(set); break;
    case MODE : Change_General_Mode(set); break;
    case ACCU : Change_Accuracy(set); break;
    case PLTP : Change_Plot_Parameters(set); break;
    case ROTA : Change_Contour_Rotation(set); break;
    case VOLG : Change_Volume_Grid(set); break;
    case BNDD : Change_Bonding(set); break;
    default   : break;
  }
}

/**************************************************************************/
/**************************************************************************/

int Change_Origin_Atom(Mol *M, Set *set, char *line)
{
  int new=0,old=Get_Atom_Number(M->origin);
  Atms *NewAt=NULL;

  if(M==NULL) return(0);
  if(!strlen(line))
  {
    sprintf(line,"Enter new origin atom number or label (%d|%s):",old,M->origin->name);
    Get_Input_Line(line,set->input);
  }
  if(strlen(line))
  {
    if(((NewAt=Goto_Atom_Name(M->atoms,line))==NULL)
     &&((new=Int_String(line))==0)) Error_Message(E_UNKNOW,"Change Origin Atom");
    else
    {
      if(NewAt) M->origin=NewAt;
      else if(new) M->origin=Goto_Atom(M->atoms,new);
      if(M->origin)
      {
        Calculate_Parameters(M,M->origin->v,set);
        End_Fractional_Coordinates(M);
      }
    }
  }
  return(1);
}

/**************************************************************************/

int Change_Origin_Group(Mol *M, Set *set, char *line)
{
  int new=0,old=0,num;
  Grps *group=NULL;
  char *args[MAXARG];
  char inline[LINELEN];
  char outline[LINELEN];

  if(M==NULL) return(0);
  if(M->groups==NULL)
  {
    Error_Message(E_NOGRPS,"Change Origin Group");
    return(0);
  }
  if(M->origin!=NULL)
  {
    old=M->origin->group;
    if((group=Goto_Group(M->groups,old))==NULL)
    {
      Error_Message(E_NOFGRP,"Change Origin Group");
      return(0);
    }
  }
  if(strlen(line)<1)
  {
    sprintf(line,"Enter new origin group number or label (%d|%s):",old,group->name);
    Get_Input_Line(line,set->input);
  }
  if((strlen(line))
   &&(strcpy(inline,line)!=NULL)
   &&((num=Get_Group_Names(inline,outline,args))>0))
  {
    if(((group=Goto_Group_Name(M->groups,args[0]))==NULL)
     &&((new=Int_String(args[0]))==0)) Error_Message(E_UNKGRP,"Change Origin Group");
    else
    {
      if((group!=NULL)&&(group->main!=NULL)) M->origin=group->main;
      else if((new)&&((group=Goto_Group(M->groups,new))!=NULL)&&(group->main!=NULL))
        M->origin=group->main;
      if(M->origin)
      {
        Calculate_Parameters(M,M->origin->v,set);
        End_Fractional_Coordinates(M);
      }
    }
  }
  return(1);
}

/**************************************************************************/

int Change_Origin_Distance(Mol *M, Set *set, char *line)
{
  double distance=0.0;
  struct vector V={0.0,0.0,0.0};

  if(M==NULL) return(0);
  if(!strlen(line))
  {
    sprintf(line,"Enter new origin distance from bonded atom (%s: number %d) (%f):"
                ,M->main_atom->name,Get_Atom_Number(M->main_atom)
                ,M->main_atom->distance);
    Get_Input_Line(line,set->input);
  }
  if (strlen(line))
  {
    sscanf(line,"%lf",&distance);
    distance-=M->main_atom->distance;
    V=SxV(-distance,unit_vector(M->main_atom->v));
    Calculate_Parameters(M,V,set);
    End_Fractional_Coordinates(M);
  }
  return(1);
}

/**************************************************************************/

void Change_Molecule_Origin(char *args, char *line, Mol *M, Set *set)
{
  if(M!=NULL)
  {
    switch(args[0])
    {
      case ATOM : Change_Origin_Atom(M,set,line); break;
      case GRPO : Change_Origin_Group(M,set,line); break;
      case DIST : Change_Origin_Distance(M,set,line); break;
      default   : break;
    }
  }
}

/**************************************************************************/
/**************************************************************************/

int Change_Molecule_Group(Mol *M, Set *set, char *line)
{
  Grps *group=NULL;
  char *args[MAXARG];
  int num,i;
  char newline[LINELEN];
  char name[LINELEN];
  char atex[LINELEN];
  double volrad=0.0;

  if(M==NULL) return(0);
  if(M->groups==NULL)
  {
    Error_Message(E_NOGRPS,"Change Molecule Groups");
    return(0);
  }
  if((num=Get_Arguements(line,args))<1)
  {
    sprintf(newline,"Enter name of group to modify :");
    Get_Input_Line(newline,set->input);
    sscanf(newline,"%s",name);
  }
  else sscanf(args[0],"%s",name);
  name[GRP_LEN]=0;
  
  if((group=Goto_Group_Name(M->groups,name))==NULL)
  {
    Error_Message(E_UNKGRP,"Change Molecule Group");
    return(0);
  }

  if(num<2)
  {
    sprintf(newline,"Enter new radius for cavity volume [%f]:",group->volrad);
    Get_Input_Line(newline,set->input);
    sscanf(newline,"%lf",&volrad);
  }
  else sscanf(args[1],"%lf",&volrad);
  if(volrad<0) volrad=0;
  if(volrad>MAXVRAD) volrad=MAXVRAD;

  atex[0]=0;
  if(num>2) for(i=2;iatex);

  for(;group!=NULL;group=group->next)
  {
    if(strcmp(group->name,name)!=0) continue;
    sprintf(line,"Changing group %s -> volrad = %f aliases = %s",name,volrad,atex);
    Out_Message(line,O_NEWLN);
    group->volrad=volrad;
    strcpy(group->atex,atex);
  }

  return(1);
}

/**************************************************************************/

int Change_Molecule_Plane(Mol *M, Set *set, char *line)
{
  double val[10];
  int num=0,i=0;
  Vector Y;

  if(M==NULL) return(0);
  for(i=0;i<10;i++) val[i]=0.0;
  if(!strlen(line))
  {
    sprintf(line,"Enter projection plane definition (0 to 10 doubles):");
    Get_Input_Line(line,set->input);
  }
  if(strlen(line))
  {
    num=sscanf(line,"%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
      &val[0],&val[1],&val[2],&val[3],&val[4],&val[5],&val[6],&val[7],&val[8],&val[9]);
    for(i=0;i100.0))
      {
        sprintf(line,"val = %f in Change Molecule Plane",val[i]);
        Error_Message(E_OUTRNG,line);
      }
      else switch(i)
      {
        case 9 : M->plane_maxP=val[i]; break;
        case 8 : M->plane_minP=val[i]; break;
        case 7 : M->plane_maxT=val[i]; break;
        case 6 : M->plane_minT=val[i]; break;
        case 5 : M->plane_x.z=val[i]; break;
        case 4 : M->plane_x.y=val[i]; break;
        case 3 : M->plane_x.x=val[i]; break;
        case 2 : M->plane.z=val[i]; break;
        case 1 : M->plane.y=val[i]; break;
        case 0 : M->plane.x=val[i]; break;
        default: break;
      }
    }
  }
  if(AlmostZero(Vmag(M->plane))) M->plane=Vequal(0.0,0.0,1.0);
  if(AlmostZero(Vmag(M->plane_x))) M->plane=Vequal(1.0,0.0,0.0);
  M->plane=unit_vector(M->plane);
  M->plane_x=unit_vector(M->plane_x);
  num=0;
  while((num<3)&&(AlmostZero(VangleV(M->plane,M->plane_x))))
  {
    Rotate_Indices_Right(&M->plane_x);
    num++;
  }
  if(num==3)
  {
    M->plane=Vequal(0.0,0.0,1.0);
    M->plane_x=Vequal(1.0,0.0,0.0);
  }
  M->plane_x=unit_vector(M->plane_x);
  Y=unit_vector(cross_product(M->plane,M->plane_x));
  M->plane_x=cross_product(Y,M->plane);
  sprintf(line,"V=(%4.2f,%4.2f,%4.2f) X=(%4.2f,%4.2f,%4.2f)",
    M->plane.x,M->plane.y,M->plane.z,
    M->plane_x.x,M->plane_x.y,M->plane_x.z);
  Out_Message(line,O_NEWLN);
  sprintf(line,"Theta=[%5.3f;%5.3f] Phi=[%5.3f;%5.3f]",
    M->plane_minT,M->plane_maxT,
    M->plane_minP,M->plane_maxP);
  Out_Message(line,O_NEWLN);
  return(1);
}

/**************************************************************************/

void Change_Molecule_Data(char *args, Mol *M, Set *set, char *line)
{
  if(M!=NULL)
  {
    switch(args[0])
    {
      case CMGP : Change_Molecule_Group(M,set,line);
      case PLAN : Change_Molecule_Plane(M,set,line);
      default   : break;
    }
  }
}

/**************************************************************************/
/**************************************************************************/

void Change_Parameter_Table(char *args, Set *set)
{
  Error_Message(E_NOTAVL,"Change Parameter Table");
}

/**************************************************************************/
/********************** view data, settings and parameters ****************/
/**************************************************************************/

void View_General_Settings(Set *set)
{
  char line[100];

  sprintf(line,"General profile array size = %d",set->size);
  Out_Message(line,O_NEWLN);
  sprintf(line,"Area versus theta profile array size = %d",set->t_size);
  Out_Message(line,O_NEWLN);
  sprintf(line,"Area versus phi profile array size = %d",set->p_size);
  Out_Message(line,O_NEWLN);
  sprintf(line,"Numerical solid angle array size = %d",set->n_size);
  Out_Message(line,O_NEWLN);
  sprintf(line,"Radial profile range from %5.2f to %5.2f Anstroms",set->min,set->max);
  Out_Message(line,O_NEWLN);
  sprintf(line,"Angular profile range from %5.2f to %5.2f Anstroms",set->a_min,set->a_max);
  Out_Message(line,O_NEWLN);
  sprintf(line,"Minimum error for Simpsons rule integration is %f",set->eps);
  Out_Message(line,O_NEWLN);
  if(set->mode&RSET_BIT)
    sprintf(line,"Above general ranges used intead of initial molecule ranges");
  else if(set->mode&RDEF_BIT)
    sprintf(line,"Initial molecule ranges used intead of above general ranges used");
  else sprintf(line,"Previous molecule ranges used intead of above general ranges used");
  Out_Message(line,O_NEWLN);
  if(set->mode&MOVG_BIT) sprintf(line,"New G positioning used in old leach solid angle"); Out_Message(line,O_NEWLN);
  if(set->mode&SCOR_BIT) sprintf(line,"Faulty segment correction used in old leach solid angle"); Out_Message(line,O_NEWLN);
  if(set->mode&NBO_BIT)  sprintf(line,"Only non-bonded overlaps used in VAO and SAO"); Out_Message(line,O_NEWLN);
  if(set->mode&NGR_BIT)  sprintf(line,"Only inter-group overlaps used in VAO and SAO"); Out_Message(line,O_NEWLN);
  if(set->mode&CONT_BIT) sprintf(line,"Angular-radial contour data output in numerical solid angle"); Out_Message(line,O_NEWLN);
  if(set->mode&PERS_BIT) sprintf(line,"Perspective mode used in contour plots"); Out_Message(line,O_NEWLN);
  if(set->mode&VDW_BIT)  sprintf(line,"Van der Waals radii used for all atoms"); Out_Message(line,O_NEWLN);
  if(set->mode&MTC_BIT)  sprintf(line,"Monte Carlo approach used in numerical volume"); Out_Message(line,O_NEWLN);

  if(set->modev&UCAV_BIT) sprintf(line,"Dummy unit cell atoms are visible"); Out_Message(line,O_NEWLN);
  if(set->modev&ABOX_BIT) sprintf(line,"Unit cell automatically boxed in cavity volume"); Out_Message(line,O_NEWLN);
  if(set->modev&REDU_BIT) sprintf(line,"Redundant atoms are automatically removed after box"); Out_Message(line,O_NEWLN);
  if(set->modev&SNAM_BIT) sprintf(line,"Atom names are used in duplication test"); Out_Message(line,O_NEWLN);
  if(set->modev&TWOD_BIT) sprintf(line,"Atom listing includes 2D coordinates"); Out_Message(line,O_NEWLN);
  if(set->modev&FRAX_BIT) sprintf(line,"Atom listing includes fractional coordinates"); Out_Message(line,O_NEWLN);
  if(set->modev&POLA_BIT) sprintf(line,"Atom listing includes polar coordinates"); Out_Message(line,O_NEWLN);
  if(set->modev&FADE_BIT) sprintf(line,"Unused atoms faded in orthogonal plot"); Out_Message(line,O_NEWLN);
}

/**************************************************************************/

void View_Parameter_Table(Parm *parm)
{
  char line[LINELEN];
  int n;

  sprintf(line,"Table of atomic parameters");
  Out_Message(line,O_NEWLN);
  for(parm=First_Parm(parm);parm!=NULL;parm=parm->next)
  {
    sprintf(line,"%-3s : v.d.Waals radius = %f  Covalent radius = %f"
                ,parm->name,parm->Vrad,parm->Crad);
    Out_Message(line,O_BLANK);
    for(n=0;nalias[n][0]!=0)
      {
        sprintf(line,", alias %s",parm->alias[n]);
        Out_Message(line,O_BLANK);
      }
    }
    Out_Message("",O_NEWLN);
  }
}

/**************************************************************************/

void View_Value(Ster *ster, double value, char *string, char *units)
{
  char line[LINELEN];
  if ((ster!=NULL)&&(value!=0.0))
  {
    sprintf(line,"%s%s = %f %s"
          ,string,ster->name,value,units);
    Out_Message(line,O_NEWLN);
  }
}

/**************************************************************************/

void View_Cone_Value(Ster *ster, double value, char *string)
{
  char line[LINELEN];
  if ((ster!=NULL)&&(value!=0.0))
  {
    sprintf(line,"%s %s = %f radians, %f degrees, %6.3f circle"
          ,string,ster->name,value,RtoD(value)
          ,value/(2.0*PI));
    Out_Message(line,O_NEWLN);
  }
}

/**************************************************************************/

void View_Solid_Value(Ster *ster, double value, char *string)
{
  char line[LINELEN];
  if ((ster!=NULL)&&(value!=0.0))
  {
    sprintf(line,"%s %s = %f steradians, %f degrees, %6.3f sphere"
          ,string,ster->name,value,StoD(value)
          ,value/(4.0*PI));
    Out_Message(line,O_NEWLN);
  }
}

/**************************************************************************/

int View_Molecule_Symmetry(Mol *M)
{
  char line[LINELEN];
  char tstr[10];
  Symm *s=NULL;

  if(M==NULL) return(0);
  if(M->symmetry==NULL) return(0);
  sprintf(line,"Symmetry and centring operators for current molecule");
  Out_Message(line,O_NEWLN);
  for(s=First_Symm(M->symmetry);s!=NULL;s=s->next)
  {
    if(s->mode&CENT_BIT) strcpy(tstr,"CENT"); else strcpy(tstr,"SYMM");
    sprintf(line,"%4s %6.2f%6.2f%6.2f%8.3f   %6.2f%6.2f%6.2f%8.3f   %6.2f%6.2f%6.2f%8.3f",tstr
                ,s->M.x.x,s->M.x.y,s->M.x.z,s->t.x
                ,s->M.y.x,s->M.y.y,s->M.y.z,s->t.y
                ,s->M.z.x,s->M.z.y,s->M.z.z,s->t.z);
    Out_Message(line,O_NEWLN);
  }
  return(1);
}

/**************************************************************************/

void View_Molecule_Atoms(Mol *M, Set *set)
{
  Atms *at=NULL;
  int count=0;
  char mark[2];
  char imp[2];
  char line[LINELEN];
  
  sprintf(line,"Atomic Data:");
  Out_Message(line,O_NEWLN);
  Out_Message("No. type:Name  radius (   x  ,   y  ,   z  )",O_BLANK);
  if(set->modev&TWOD_BIT)
    Out_Message(" (   X  ,   Y  )",O_BLANK);
  if(set->modev&FRAX_BIT)
    Out_Message(" (  fx  ,  fy  ,  fz  )",O_BLANK);
  if(set->modev&POLA_BIT)
    Out_Message(" (   d  , theta,  phi )",O_BLANK);
  Out_Message(" SVangle",O_NEWLN);
  Out_Message("==========================================================================",O_NEWLN);

  count=0;

  for(at=First_Atom(M->atoms);at!=NULL;at=at->next)
  {
    count++; if(count>M->num_atoms) Error_Message(E_ATMNUM,"View Molecule Atoms");
    if((!(set->modev&UCAV_BIT))&&(at->group==FULL_BIT)) continue;
    mark[0]=' ';
    mark[1]=' ';
    imp[0]=' ';
    imp[1]=' ';
    if(at==M->origin) imp[0]='o';
    if(count==Get_Atom_Number(M->main_atom)) imp[1]='*';
    if(!(at->stat&MAIN_BIT)) { mark[0]='<'; mark[1]='>'; }
    sprintf(line,"%2d:%d%c%c%c%3s:%-4s %5.3f (%6.3f,%6.3f,%6.3f)"
		,count,at->group,imp[0],imp[1],mark[0],at->type,at->name,at->radius
		,at->v.x,at->v.y,at->v.z);
    Out_Message(line,O_BLANK);
    if(set->modev&TWOD_BIT)
    {
      sprintf(line," (%6.3f,%6.3f)",at->tv.x,at->tv.y);
      Out_Message(line,O_BLANK);
    }
    if(set->modev&FRAX_BIT)
    {
      sprintf(line," (%6.3f,%6.3f,%6.3f)",at->fv.x,at->fv.y,at->fv.z);
      Out_Message(line,O_BLANK);
    }
    if(set->modev&POLA_BIT)
    {
      sprintf(line," (%6.3f,%6.3f,%6.3f)",at->distance,at->theta,at->phi);
      Out_Message(line,O_BLANK);
    }
    sprintf(line," %5.3f %c ",at->SVangle,mark[1]);
    Out_Message(line,O_NEWLN);
  }
}

/**************************************************************************/

void View_Molecule_Groups(Mol *M)
{
  int count=0;
  char line[LINELEN];
  Grps *group=NULL;
  
  for(count=1,group=First_Group(M->groups);group!=NULL;group=group->next,count++)
  {
    sprintf(line,"Group %d:%s (atom %d:%s of %d)"
                ,count,group->name,Get_Atom_Number(group->main)
                ,group->main->name,group->num_atoms);
    if(!AlmostZero(group->SVangle))
    {
      Out_Message(line,O_BLANK);
      sprintf(line," ; max SV angle = %6.3f radians",group->SVangle);
    }
    if(!AlmostZero(group->volume))
    {
      Out_Message(line,O_BLANK);
      sprintf(line," ; group volume =%8.2f cubic angstroms",group->volume);
    }
    if(!AlmostZero(group->cavity))
    {
      Out_Message(line,O_BLANK);
      sprintf(line," ; cavity volume =%8.2f cubic angstroms",group->cavity);
    }
    if(!AlmostZero(group->radvol))
    {
      Out_Message(line,O_BLANK);
      sprintf(line," ; radial volume =%8.2f cubic angstroms",group->radvol);
    }
    if(!AlmostZero(group->radcav))
    {
      Out_Message(line,O_BLANK);
      sprintf(line," ; radcav volume =%8.2f cubic angstroms",group->radcav);
    }
    if(!AlmostZero(group->volrad))
    {
      Out_Message(line,O_BLANK);
      sprintf(line," ; cavity radius =%6.2f angstroms",group->volrad);
    }
    Out_Message(line,O_NEWLN);
  }
}

/**************************************************************************/

int View_Molecule_Name(Mol *M)
{
  char line[LINELEN];
  int n,a,g,s;
  Atms *A=NULL;
  Grps *G=NULL;
  Ster *S=NULL;
  
  if(M==NULL) return(0);
  n=Get_Molecule_Number(M);
  for(a=0,A=First_Atom(M->atoms)  ;A!=NULL;A=A->next,a++);
  for(g=0,G=First_Group(M->groups);G!=NULL;G=G->next,g++);
  for(s=0,S=First_Steric(M->ster) ;S!=NULL;S=S->next,s++);
  sprintf(line,"Molecule %-*s is number %3d with %4d atoms in %2d groups, having %1d steric results"
              ,10,M->name,n,a,g,s);
  Out_Message(line,O_NEWLN);
  return(1);
}

/**************************************************************************/

int View_Molecule_Data(Mol *M, char *args, Set *set)
{
  int count=0;
  char perc='%';
  char line[LINELEN];
  Ster *ster=NULL;
  
  if(M==NULL) return(0);
  sprintf(line,"Data for molecule %s.",M->name);
  Out_Message(line,O_NEWLN);
  sprintf(line,"Number of atoms = %d",M->num_atoms);
  Out_Message(line,O_NEWLN);
  sprintf(line,"minR=%f maxR=%f",M->minR,M->maxR);
  Out_Message(line,O_NEWLN);

  if(M->mode&FRAC_BIT)
  {
    sprintf(line,"Unit Cell Parameters: a=%7.3f b=%7.3f c=%7.3f alpha=%8.3f beta=%8.3f gamma=%8.3f"
                ,M->a,M->b,M->c,M->al,M->be,M->ga);
    Out_Message(line,O_NEWLN);
    sprintf(line,"Unit Cell Volume = %f cubic angstroms",M->cellvol);
    Out_Message(line,O_NEWLN);
    sprintf(line,"Number of general equivalent positions is %d",M->geneq);
    Out_Message(line,O_NEWLN);
  }

  sprintf(line,"Original Tolman cone angle = %f radians, %f degrees, %3.0f%c circle"
              ,M->tolman,RtoD(M->tolman)
              ,100*M->tolman/(2.0*PI),perc);
  Out_Message(line,O_NEWLN);

  for(ster=First_Steric(M->ster);ster!=NULL;ster=ster->next)
  {
    if(ster->type==S_VOLM)
      View_Value(ster,ster->tot_val,"","cubic angstroms");
    else if(ster->type==S_PROJ)
      View_Value(ster,ster->tot_val,"","square angstroms");
    else if((ster->type==S_ANGU)||(ster->type==S_CONE)
     ||(ster->type==S_TOLM)||(ster->type==S_VAOV))
    {
      View_Cone_Value(ster,ster->tot_val,"Total");
      View_Cone_Value(ster,ster->max_val,"Maximum profile");
      View_Value(ster,ster->pr_area,"Profile Area ","radian.angstroms");
      View_Cone_Value(ster,ster->tot_con,"Conformer averaged");
    }
    else
    {
      View_Solid_Value(ster,ster->tot_val,"Total");
      View_Solid_Value(ster,ster->max_val,"Maximum profile");
      View_Value(ster,ster->pr_area,"Profile Area ","steradian.angstroms");
      View_Solid_Value(ster,ster->tot_con,"Conformer averaged");
    }
    View_Value(ster,ster->peak_R,"Radius at profile peak for ","angstroms");
  }

  for(count=0;args[count]!=0;count++)
  {
    switch(args[count])
    {
      case ATMV : View_Molecule_Atoms(M,set); break;
      case GRPV : View_Molecule_Groups(M); break;
      case SYMV : View_Molecule_Symmetry(M); break;
      default   : break;
    }
  }

  return(1);
}

/**************************************************************************/
/**********************   other functions  ********************************/
/**************************************************************************/

void Initialize_Settings(Set *set)
{
  set->parm=NULL;          /* pointer to table of parameters              */
  set->input=NULL;         /* text input stream                           */
  set->contour=NULL;       /* contou output file                          */
  set->rotation=DEF_ROT;   /* contour plot rotation angle about z axis    */
  set->multi=MAX_OVER;     /* number of overlaps to consider in craig SA  */
  set->size=P_DEF;         /* size of general profile array               */
  set->t_size=AT_DEF;      /* area vs. theta profile size                 */
  set->p_size=AT_DEF;      /* area vs. phi profile size                   */
  set->n_size=S_DEF;       /* size of array used for numerical solid angle*/
  set->c_size=C_DEF;       /* size of array used for graphical 3D profile */
  set->vx=V_DEF;           /* size of x array for numerical volume        */
  set->vy=V_DEF;           /* size of y array for numerical volume        */
  set->vz=V_DEF;           /* size of z array for numerical volume        */
  set->min=100.0;
  set->max=0.0;            /* minimum and maximum radial profile range    */
  set->a_min=0.0;
  set->a_max=2*PI;         /* minimum and maximum angular profile range   */
  set->eps=E_PS;           /* minimum error in simpsons rule              */
  set->veps=VE_PS;         /* minimum error in MC volume calculation      */
  set->tgap=VT_DEF;        /* default time gap between MCV outputs        */
  set->avmmcv=AVM_MCV;     /* point at which to start averaging in MCV    */
  set->minmcv=MIN_MCV;     /* minimum count for Monte Carlo Volume        */
  set->maxmcv=MAX_MCV;     /* max grid points to test in MCV              */
  set->smooth=SM_DEF;      /* profile smoothing window                    */
  set->vrcut=VRC_DEF;      /* volrad cutoff in profile                    */
  set->bond_minf=BOND_MINF;/* fractional min. bond length in find bond    */
  set->bond_maxf=BOND_MAXF;/* fractional max. bond length in find bond    */
  set->mode=RSET_BIT       /* use general settings profile range          */
           |PERS_BIT       /* use perspective mode in contours            */
           |MTC_BIT;       /* use monte carlo volume calculation          */
  set->modev=0;            /* general viewing mode                        */
  set->plot=NULL;          /* pointer to plot pipe to plotting program    */
  set->x_ran=X_RAN;
  set->y_ran=Y_RAN;        /* x and y ranges for perspective plot         */
}

/**************************************************************************/

Mol *Read_Arguements(int argc, char *argv[], Mol *M, Set *set)
{
  int count;

  for(count=1;count");
  else sprintf(newline,"%d:Steric<%s>",Get_Molecule_Number(M),M->name);
  if((n=Get_Command_Line(newline,set->input))==0)
  {
    fclose(set->input);
    set->input=stdin;
  }
  if(n==2) return(EXIT);
  if(newline[0]=='#') return(BLANK);
  if(newline[0]=='>')
  {
    if(output!=stdout) fclose(output);
    sscanf(newline+1,"%s",line);
    if(!strlen(line))
    {
      output=stdout;
      Out_Message("Output redirected to stdout",O_NEWLN);
    }
    else if((output=fopen(line,"wt"))==NULL)
    {
      output=stdout;
      Error_Message(E_OPNOUT,"Get Command");
    }
    else printf("Output redirected to file %s\n",line);
    line[0]=0;
    return(0);
  }
  if(newline[0]=='!')
  {
    strcpy(line,newline+1);
    return(0);
  }
  if(newline[0]=='&')
  {
    set->mode=set->mode|ALL_BIT;
    i=1;
  }
  else
  {
    set->mode=set->mode&(FULL_BIT^ALL_BIT);
    i=0;
  }
  if(strchr(newline,'\n')!=NULL) *strchr(newline,'\n')=0;
  strcpy(line,newline+i);
  if(Find_Commands(args,line,newline+i,comm)) return(args[0]);
  else return(0);
}

/**************************************************************************/

void Free_Memory(Comm *comm, Parm *parm, Mol *M)
{
  if(comm!=NULL) Close_All_Commands(comm);
  if(parm!=NULL) free(parm);
  M=First_Molecule(M);
  while(M!=NULL) M=Close_Current_Molecule(M);
}

/**************************************************************************/
/********************** main menu operations ******************************/
/**************************************************************************/

Mol *File_Operations(char *args, Mol *molecule, char *line, Set *set)
{
  Mol *M=molecule;
  switch(args[0])
  {
    case LOAD : molecule=Load_Molecules(line,molecule,set); break;
    case SAVE : if(!set->mode&ALL_BIT) Save_Molecules(args,molecule,line);
                else for(M=First_Molecule(molecule);M!=NULL;M=M->next)
                {
                  Save_Molecules(args,M,line);
                }
                break;
    case SVAS : Save_Molecules(args,molecule,line); break;
    case SHEL : Open_Shell(); break;
    default   : break;
  }
  return(molecule);
}

/**************************************************************************/

Mol *Molecule_Operations(char *args, Mol *molecule, char *line, Set *set)
{
  Mol *M=molecule;
  Mol *old=NULL;
  if(molecule==NULL) return(NULL);
  switch(args[0])
  {
    case VIEW : if(!(set->mode&ALL_BIT)) View_Molecule_Data(molecule,"",set);
                else for(M=First_Molecule(molecule);M!=NULL;M=M->next)
                {
                  View_Molecule_Data(M,"",set);
                }
                break;
    case NEXT : molecule=Next_Molecule(molecule); break;
    case PREV : molecule=Previous_Molecule(molecule); break;
    case GOTO : molecule=Goto_Particular_Molecule(molecule,line,set); break;
    case KILL : if(!(set->mode&ALL_BIT)) Kill_Molecule_Data(args+1,line,molecule,set);
                else for(M=First_Molecule(molecule);M!=NULL;M=M->next)
                {
                  Kill_Molecule_Data(args+1,line,M,set);
                }
                break;
    case CLOS : if(!(set->mode&ALL_BIT)) molecule=Test_Close_Current_Molecule(molecule,line);
                else for(M=First_Molecule(molecule);M!=NULL;)
                {
                  old=M;
                  M=Test_Close_Current_Molecule(M,line);
                  if(old==M) M=M->next;
                  if(M!=NULL) molecule=M;
                }
                break;
    default   : break;
  }
  return(molecule);
}

/**************************************************************************/

void Symmetry_Operations(char *args, Mol *molecule, char *line, Set *set)
{
  Mol *M=NULL;
  if(molecule!=NULL)
  {
    for(M=First_Molecule(molecule);M!=NULL;M=M->next)
    {
      if(set->mode&ALL_BIT) molecule=M;
      switch(args[0])
      {
        case TRAN : Transform_Data(molecule,line,set); break;
        case XPND : Expand_Data(molecule,set); break;
        case BOXU : Box_Atoms(molecule,line,set); break;
        case BOND : Find_All_Bonds(molecule,set,0,VIS_BIT); break;
        case GRUP : Find_All_Bonded_Groups(molecule,0); break;
        case REDU : Remove_Duplicate_Atoms(molecule,set); break;
        case XGRR : Exclude_Atoms_Outside_Groups(molecule,set,line); break;
        case SHLV : Calculate_Group_Shell(args+1,molecule,set,line); break;
        case SVOL : Calculate_Crystal_Volumes(args+1,molecule,set,line); break;
        case SPRO : Calculate_Crystal_Areas(args+1,molecule,set,line); break;
        default   : break;
      }
      if(!(set->mode&ALL_BIT)) break;
    }
  }
}

/**************************************************************************/

void Calculations(char *args, Mol *molecule, Set *set, char *line)
{
  Mol *M=NULL;
  char fname[LINELEN];
  if(molecule!=NULL)
  {
    for(M=First_Molecule(molecule);M!=NULL;M=M->next)
    {
      if(set->mode&ALL_BIT) molecule=M;
      if(set->mode&CONT_BIT)
      {
        strcpy(fname,molecule->Fname);
        if(fname[0]==0) sprintf(fname,"M%d",Get_Molecule_Number(molecule));
        New_Extension(fname,".ctr");
        if(set->contour!=NULL) fclose(set->contour);
        set->contour=fopen(fname,"wt");
      }
      switch(args[0])
      {
        case TOTA : Calculate_Steric_Value(args+1,molecule,set); break;
        case CONF : Calculate_Steric_Conformer_Average(args+1,molecule,set); break;
        case PROF : Calculate_Steric_Profile(args+1,molecule,set); break;
        case VOLM : Calculate_Volumes(args,molecule,set,line); break;
        case AREA : Calculate_Areas(args+1,molecule,set); break;
        default   : break;
      }
      if(set->contour!=NULL) { fclose(set->contour); set->contour=NULL; }
      if(!(set->mode&ALL_BIT)) break;
    }
  }
}

/**************************************************************************/

int Plot(char *args, Mol *molecule, Set *set, char *line)
{ 
  Mol *M=NULL;
  char ext[10];
  char fname[100];
  int ans=0;
  if(molecule!=NULL)
  {
    for(M=First_Molecule(molecule);M!=NULL;M=M->next)
    {
      if(set->mode&ALL_BIT) molecule=M;
      switch(args[0])
      {
        case ORTH : Draw_Projected_View(molecule,set); break;
        case MOLP : Draw_Special_View(molecule,set); break;
        case CART : ans=Draw_Steric_Profile(args+1,molecule,set,args[0]); break;
        case POLA : ans=Draw_Steric_Profile(args+1,molecule,set,args[0]); break;
        case STER : Draw_Steric_Parameter(args+1,molecule,set); break;
        case CONT : Draw_Molecule_Contour(molecule,set,CONT_SCRIPT); break;
        case PRFC : Draw_Molecule_Contour(molecule,set,PROF_SCRIPT); break;
        default   : if(set->plot!=NULL) fprintf(set->plot,"%s\n",line); break;
      }
      if(!ans) break;
      if(!(set->mode&ALL_BIT)) break;
    }
    if(ans)
    {
      if(args[0]==POLA)
      {
        fprintf(set->plot,"set polar\n");
        strcpy(ext,".ang");
      }
      else
      {
        fprintf(set->plot,"set nopolar\n");
        strcpy(ext,".rad");
      }
      fprintf(set->plot,"plot");
      for(M=First_Molecule(molecule);M!=NULL;M=M->next)
      {
        if(set->mode&ALL_BIT) molecule=M;
        strcpy(fname,molecule->Fname);
        New_Extension(fname,ext);
        fprintf(set->plot," \"%s\" with lines",fname); 
        if(!(set->mode&ALL_BIT)) break;
        else if(M->next!=NULL) fprintf(set->plot,", ");
      }
      fprintf(set->plot,"\n");
      fflush(set->plot);
    }
  }
  return(ans);
}

/**************************************************************************/

void Change(char *args, char *line, Mol *molecule, Set *set)
{
  Mol *M=NULL;
  switch(args[0])
  {
    case SETT : Change_General_Settings(args+1,set); break;
    case ORIG : if(!(set->mode&ALL_BIT)) Change_Molecule_Origin(args+1,line,molecule,set);
                else for(M=First_Molecule(molecule);M!=NULL;M=M->next)
                {
                  Change_Molecule_Origin(args+1,line,M,set);
                }
                 break;
    case MOLE : if(!(set->mode&ALL_BIT)) Change_Molecule_Data(args+1,molecule,set,line);
                else for(M=First_Molecule(molecule);M!=NULL;M=M->next)
                {
                  Change_Molecule_Data(args+1,M,set,line);
                } break;
    case PARA : Change_Parameter_Table(args+1,set); break;
    default   : break;
  }
}

/**************************************************************************/

void View(char *args, Mol *molecule, Set *set)
{
  Mol *M=NULL;
  switch(args[0])
  {
    case SETT : View_General_Settings(set); break;
    case MOLE : if(molecule==NULL) break;
                else if(!(set->mode&ALL_BIT)) View_Molecule_Data(molecule,args+1,set);
                else for(M=First_Molecule(molecule);M!=NULL;M=M->next)
                {
                  View_Molecule_Data(M,args+1,set);
                }
                break;
    case PARA : View_Parameter_Table(set->parm); break;
    default   : for(M=First_Molecule(molecule);M!=NULL;M=M->next)
                  View_Molecule_Name(M);
                break;
  }
}

/**************************************************************************/

int Help(char *argline, Comm *comm)
{
  char args[LINELEN];
  char line[LINELEN];
  char fname[LINELEN];
  FILE *hlp=NULL;
  char *argptr[MAXARG];
  int num=0,count=0,i;

  Find_Commands(args,line,argline,comm);
  args[MAXARG-1]=0;
  num=strlen(args);

  sprintf(fname,"steric.hlp");
  if((hlp = fopen(fname, "rt")) == NULL)
  {
    sprintf(fname,"%s/steric.hlp",STERICHOME);
    if((hlp = fopen(fname, "rt")) == NULL)
    {
      sprintf(fname,"Help: [steric.hlp] | [%s/steric.hlp]",STERICHOME);
      Error_Message(E_NOHELP,fname);
      return(0);
    }
  }
  if(num>0)
  {
    strcpy(line,"Help: ( ");
    for(count=strlen(line),num+=count,i=0;count+i0)
      for(count=0;countplot!=NULL) fclose(set->plot);
  Free_Memory(comm,set->parm,molecule);
  Error_Message(error,"");
}

/**************************************************************************/
/******************  Main program controls user  **************************/
/******************  interface and all function  **************************/
/******************            calls             **************************/
/**************************************************************************/

int main(int argc, char *argv[])
{
  Set set;                  /* structure containing general settings      */
  Mol *molecule=NULL;       /* pointer to current data set                */
  Comm *first_command=NULL; /* pointer to branched list of commands       */
  char line[LINELEN];       /* pointer to command line entered            */
  char args[LINELEN];       /* pointer to arguements of command line      */
  char choice=0;            /* single byte command                        */
  char name[30];

  output=stdout;
  Title();
  Initialize_Settings(&set);
  if((first_command=Initialize_Commands())==NULL)
    Soft_Exit(first_command,&set,molecule,E_NOCOMM);
  Load_Parameters(&set,(FILE *)NULL,(char *)NULL);
  if(argc>1) molecule=Read_Arguements(argc,argv,molecule,&set);
  if(set.input==NULL)
  {
    if((set.input=fopen("steric.ini","rt"))==NULL)
    {
      sprintf(line,"%s/steric.ini",STERICHOME);
      if((set.input=fopen(line,"rt"))==NULL) set.input=stdin;
    }
  }
  do
  {
    if(molecule) strcpy(name,molecule->name);
    else name[0]=0;
    switch(choice=Get_Command(molecule,line,args,first_command,&set))
    {
      case FILS : molecule=File_Operations(args+1,molecule,line,&set); break;
      case MOLE : molecule=Molecule_Operations(args+1,molecule,line,&set); break;
      case SYMM : Symmetry_Operations(args+1,molecule,line,&set); break;
      case CALC : Calculations(args+1,molecule,&set,line); break;
      case PLOT : Plot(args+1,molecule,&set,line); break;
      case CHAN : Change(args+1,line,molecule,&set); break;
      case VIEW : View(args+1,molecule,&set); break;
      case HELP : Help(line,first_command); break;
      case EXIT : Soft_Exit(first_command,&set,molecule,E_BLANK); break;
      case BLANK: line[0]=0; break;
      default   : if(line[0]) system(line); break;
    }
  } while (choice != EXIT);
  return(1);
}

/**************************************************************************/
/******************  The End ... ******************************************/
/**************************************************************************/