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

/**************************************************************************/
/******************  File manipulation functions **************************/
/**************************************************************************/
/******************        This module is        **************************/
/******************      system independant      **************************/
/**************************************************************************/

#include 
#include 
#include 
#include 

#include "sterdefn.h"      /* declaration of structures and globals       */
#include "stercomm.h"      /* definitions of all menu command options     */
#include "sterfile.h"      /* functions for all file manipulations        */
#include "crystal.h"       /* crystallographic manipulations              */
#include "stermem.h"       /* dynamic memory management module            */
#include "stertext.h"      /* all text functions for text mode            */
#include "stergrp.h"       /* group manipulation functions                */
#include "steraid.h"       /* additional functions needed                 */
#include "stercalc.h"      /* main calculation function                   */
#include "craig.h"         /* craig's multiple overlap calculation        */

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

int Initialize_Molecule(Mol *M)
{
  if(M==NULL) return(0);
  M->name[0]=0;            /* molecule name                               */
  M->Fname[0]=0;           /* molecule filename without extension         */
  M->a=1.0;                /* unit cell parameters, if appropriate        */
  M->b=1.0;
  M->c=1.0;
  M->al=90.0;
  M->be=90.0;
  M->ga=90.0;
  M->cellvol=0.0;          /* unit cell volume                            */
  M->freevol=0.0;          /* free volume in unit cell                    */
  M->geneq=1;              /* number of general equivalent positions      */
  M->centring='p';         /* type of centring condition                  */
  M->num_atoms=0;          /* number of atoms                             */
  M->multi=MAX_OVER;       /* number of overlaps in multiple overlap calc.*/
  M->main_atom=NULL;       /* atom bonded to origin defining atom (if any)*/
  M->minR=0.0;
  M->maxR=10.0;            /* min,max molecule distance from apex         */
  M->plane=Vequal(0.0,0.0,1.0);
  M->plane_x=Vequal(1.0,0.0,0.0);
  M->plane_minT=0.0;
  M->plane_maxT=PI;        /* plane rotation theta range                  */
  M->plane_minP=0.0;
  M->plane_maxP=2*PI;      /* plane rotation theta range                  */
  M->plane_T=0.0;
  M->plane_P=0.0;          /* current plane theta and phi                 */
  M->basis_x=Vequal(1.0,0.0,0.0);
  M->basis_y=Vequal(0.0,1.0,0.0);
  M->basis_z=Vequal(0.0,0.0,1.0);
  M->main_atom=NULL;       /* pointer to main atom                        */
  M->origin=NULL;          /* pointer to origin defining atom             */
  M->atoms=NULL;           /* pointer to individual atoms memory          */
  M->groups=NULL;          /* pointer to definable atomic groups          */
  M->symmetry=NULL;        /* pointer to symmetry operator memory         */
  M->ster=NULL;            /* pointer to individual calculations memory   */
  M->mode=0;               /* molecule mode                               */
  return(1);
/* next and prev pointers have already been assigned, do not fiddle !!!   */
}

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

int Initialize_Atom(Atms *atom)
{
  if(atom==NULL) return(0);
  atom->name[0]=0;         /* atom lable                                  */
  atom->type[0]=0;         /* atom type label                             */
  atom->v=Vequal(0.0,0.0,0.0);    /* position rel. to apex                */
  atom->fv=Vequal(0.0,0.0,0.0);   /* position rel. to apex                */
  atom->tv=Vequal(0.0,0.0,0.0);   /* position in 2D, on projection plane  */
  atom->distance=0.0;	   /* distance from apex                          */
  atom->theta=0.0;
  atom->phi=0.0;           /* polar position with main atom at 0,0        */
  atom->radius=0.0;        /* atomic radius                               */
  atom->SVangle=0.0;       /* effective semi-vertix angle for calc.       */
  atom->bonds=0;           /* number of bonds                             */
  atom->count=0;           /* number of times counted in craig counting   */
  atom->stat=0;  	   /* status of atom for calculations             */
  atom->group=0;           /* group numbers                               */
  atom->bond=NULL;         /* pointer to one of the atoms bonds           */
  return(1);
/* next and prev pointers have already been assigned, do not fiddle !!!   */
}

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

int Initialize_Bond(Bond *bond)
{
  if(bond==NULL) return(0);
  bond->type=0;            /* bond type (0=not bonded, 1=single, 2=double */
  bond->atom=NULL;         /* pointer to atom bonded to                   */
  return(1);
/* next and prev pointers have already been assigned, do not fiddle !!!   */
}

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

Bond *Find_Bond(Atms *A, Atms *B)
{
  Bond *bond=NULL;
  if((A==NULL)||(B==NULL)) return(NULL);
  if((bond=First_Bond(A->bond))==NULL) return(NULL);
  for(bond=First_Bond(A->bond);bond!=NULL;bond=bond->next)
    if(bond->atom==B) return(bond);
  return(NULL);
}

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

int Bond_Exists(Atms *A, Atms *B, short unsigned type)
{
  Bond *bond=NULL;
  if((A==NULL)||(B==NULL)) return(0);
  if(A->bond==NULL)return(0);
  if((bond=Find_Bond(A,B))==NULL) return(0);
  if(bond->type!=type) Error_Message(E_BDBOND,"Bond Exists");
  return(1);
}

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

void Add_New_Bond(Atms *atomA, Atms *atomB, short unsigned type)
{
  if((atomA==NULL)||(atomB==NULL)) Error_Message(E_IMBOND,"Add New Bond");
  else
  {
    if(!Bond_Exists(atomA,atomB,type))
    {
      atomA->bond=New_Bond(atomA->bond);
      Initialize_Bond(atomA->bond);
      atomA->bond->type=type;
      atomA->bond->atom=atomB;
      atomA->bonds++;
    }
    if(!Bond_Exists(atomB,atomA,type))
    {
      atomB->bond=New_Bond(atomB->bond);
      Initialize_Bond(atomB->bond);
      atomB->bond->type=type;
      atomB->bond->atom=atomA;
      atomB->bonds++;
    }
  }
}

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

void Screen_Atom_Name(char *name)
{
  int i;
  int len=strlen(name);
  for (i=0;i122)||(name[i]=='_')) name[i]=0;
  }
}

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

int Load_Parameters(Set *set, FILE *PF, char *fname)
{
  char *args[MAXARG];
  int num=0;
  Parm *new=NULL;
  int n=0;
  char name[LINELEN];
  char line[LINELEN];

  if((fname!=NULL)&&(strlen(fname)>0)) strcpy(name,fname);
  else strcpy(name,"steric.par");
  sprintf(line,"%s",name);
  if((PF==NULL)&&((PF = fopen(line,"rt")) == NULL))
  {
    sprintf(line,"%s/%s",STERICHOME,name);
    if((PF==NULL)&&((PF = fopen(line,"rt")) == NULL))
    {
      sprintf(line,"Load Parameters [%s] | [%s/%s]",name,STERICHOME,name);
      Error_Message(E_NOPARM,line);
      return (0);
    }
  }
  if(set->parm!=NULL) set->parm=Close_All_Parms(set->parm);
  while(fgets(line,LINELEN,PF)!=NULL)
  {
    if(line[0]=='#') continue;
    if(strchr(line,'\n')!=NULL) strchr(line,'\n')[0]=0;
    if(strchr(line,'#')!=NULL) strchr(line,'#')[0]=0;
    if((num=Get_Arguements(line,args))<1) continue;
    if((new=New_Parm(set->parm))==NULL) { fclose(PF); return(0); }
    set->parm=new;
    set->parm->name[0]=0;
    set->parm->Vrad=0.0;
    set->parm->Crad=0.0;
    for(n=0;nparm->alias[n][0]=0;
    for(n=0;nparm->name,name);
                  break;
        case 1  : sscanf(args[n],"%lf",&set->parm->Vrad); break;
        case 2  : sscanf(args[n],"%lf",&set->parm->Crad); break;
        default : sscanf(args[n],"%s",name);
                  name[ATM_LEN]=0;
                  strcpy(set->parm->alias[n-3],name);
                  break;
      }
    }
  }
  fclose(PF);
  set->parm=First_Parm(set->parm);
  return(1);
}

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

double Get_Atomic_Radius(int bonds, char *name, Parm *parm, unsigned mode)
{
  Parm *pa=NULL;
  int i,n;
  if(parm==NULL)
  {
    Error_Message(E_BDPARM,"Get Atomic Radius");
    return(DEF_RAD);
  }
  for(i=0;i<2;i++)
  {
    for(pa=First_Parm(parm);pa!=NULL;pa=pa->next)
    {
      if((stricmp(name,pa->name))!=0)
      {
        for(n=0;nalias[n]))==0) break;
        }
        if(n==MAX_ALS) continue;
        strcpy(name,pa->name);
      }
      if(mode&VDW_BIT) return(pa->Vrad);
      if(mode&CVL_BIT) return(pa->Crad);
      else if (bonds<2) return(pa->Vrad);
      else return(pa->Crad);
    }
    Screen_Atom_Name(name);
  }
  return(0.0);
}

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

int Calculate_Polar_Positions(Mol *M)
{
  Atms *at=NULL;
  int n=0;
  struct vector flat={0.0,0.0,0.0};
  double x,y;
  char line[LINELEN];

  if((M==NULL)||(M->atoms==NULL)) return(0);
  Set_Projected_XY(M,0.0,0.0);
  Set_Total_SVAngles(M);
  Get_Theta_Positions(M);               /* <- basis_z is calculated here */
  M->tolman=Find_Max_Group_SVangles(M);    /* M->atoms for basis_x found */
  sprintf(line,"Tolman Cone Angle is %8.4f radians",M->tolman);
  Out_Message(line,O_NEWLN);
  M->basis_x=unit_vector(cross_product(M->basis_z,unit_vector(M->atoms->v)));
  M->basis_x=unit_vector(cross_product(M->basis_x,M->basis_z));
  if(!AlmostZero(dot_product(M->basis_z,M->basis_x)))
    Error_Message(E_BDBASE,"Calculate Polar Positions");

  M->basis_y=unit_vector(cross_product(M->basis_z,M->basis_x));

/* calculate phi positions                                               */

  for(n=0,at=First_Atom(M->atoms);at!=NULL;at=at->next)
  {
    n++; if(n>M->num_atoms) Error_Message(E_ATMNUM,"Calculate Polar Position");
    flat=cross_product(M->basis_z,unit_vector(at->v));
    if(AlmostZero(at->theta)) at->phi=0.0;
    else
    {
      flat=unit_vector(flat);
      if(!AlmostZero(dot_product(flat,M->basis_z))) Error_Message(E_BDPROJ,"Calculate Polar Position");
      flat=cross_product(flat,M->basis_z);
      if(!AlmostZero(dot_product(flat,M->basis_z))) Error_Message(E_BDPROJ,"Calculate Polar Position");
      flat=unit_vector(flat);
      y=dot_product(flat,M->basis_y);
      x=dot_product(flat,M->basis_x);
      if(x>=0)
      {
        if(y>=0) at->phi=asin(y);
        else at->phi=2*PI-asin(-y);
      }
      else
      {
        if(y>=0) at->phi=PI-asin(y);
        else at->phi=PI+asin(-y);
      }
    }
  }
  return(1);
}

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

int Get_Positions(Mol *M, struct vector origin, Set *set)
{
  Atms *at=NULL;
  double distance=0.0;
  double radius=0.0;
  int n;
  double maxR=0, minR=100;

  if(M==NULL) return(0);
  if((M->mode&FRAC_BIT)&&(Vcmp(origin,Vequal(0.0,0.0,0.0))!=0))
    End_Fractional_Coordinates(M);
  for(n=0,at=First_Atom(M->atoms);at!=NULL;at=at->next)
  {
    n++; if(n>M->num_atoms) Error_Message(E_ATMNUM,"Get Positions");
    at->v=Vsum(at->v,SxV(-1,origin));
    distance=sqrt(at->v.x*at->v.x
                 +at->v.y*at->v.y
                 +at->v.z*at->v.z);
    at->distance=distance;
    radius=at->radius;
    if ((strcmp(at->name,"DUO")!=0)&&(strcmp(at->type,"DUO")!=0))
    {
      if (maxR<(distance+radius)) maxR=(distance+radius);
      if (minR>(distance-radius)) minR=(distance-radius);
    }
  }
  M->maxR=maxR;
  M->minR=minR;
  set->min=fmin(set->min,minR);
  set->max=fmax(set->max,maxR);
  return(1);
}

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

int Calculate_Parameters(Mol *M, struct vector origin, Set *set)
{
  Bond *bond=NULL;

  if((M==NULL)||(M->atoms==NULL)) return(0);
  if(M->origin) Assign_Groups(M);
  Get_Positions(M,origin,set);
  if(M->origin)
  {
    for(bond=First_Bond(M->origin->bond);bond!=NULL;bond=bond->next) 
    {
      M->main_atom=bond->atom;
      if(M->main_atom->stat&MAIN_BIT) break;
    }
    M->origin->stat=M->origin->stat&(!MAIN_BIT);
  }
  if(!(M->mode&FRAC_BIT)) Calculate_Polar_Positions(M);
  return(1);
}

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

int Find_Atomic_Radii(Mol *M, Set *set)
{
  Atms *at=NULL;
  int n=0;
  if((M==NULL)||(M->atoms==NULL)) return(0);
  at=First_Atom(M->atoms);
  for(n=0,at=First_Atom(M->atoms);at!=NULL;at=at->next)
  {
    n++; if(n>M->num_atoms) Error_Message(E_ATMNUM,"Find Atomic Radii");
    at->stat=at->stat&(FULL_BIT^GRP_BIT);
    at->radius=Get_Atomic_Radius(at->bonds,at->type,set->parm,set->mode);
    if (at->radius!=0) at->stat=MAIN_BIT;
  }
  return(1);
}

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

int Setup_Atomic_Parameters(Mol *M, Set *set)
{
  if((M==NULL)||(M->atoms==NULL)) return(0);
  Find_Atomic_Radii(M,set);
  if(M->origin) Calculate_Parameters(M,M->origin->v,set);
  else Calculate_Parameters(M,Vequal(0.0,0.0,0.0),set);
  return(1);
}

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

void Save_Steric_Atoms(FILE *MF, Mol *M)
{
  char line[LINELEN];
  int i=0,n=0;
  Atms *at=NULL;
  Bond *bond=NULL;
  for(at=First_Atom(M->atoms);at!=NULL;at=at->next)
  {
    if(at==M->origin)
    {
      i=Get_Atom_Number(at);
      strcpy(line,"ORIGIN");
    }
    else strcpy(line,"ATOM");
    fprintf(MF,"%-6s  %-*s%-*s%10.6f%10.6f%10.6f\n"
              ,line,ATM_LEN,at->name,ATM_LEN,at->type
              ,at->fv.x,at->fv.y,at->fv.z);
  }
  for(at=First_Atom(M->atoms);at!=NULL;at=at->next)
  {
    if(at->bond==NULL) continue;
    i=Get_Atom_Number(at);
    fprintf(MF,"BOND  %4d",i);
    for(bond=First_Bond(at->bond);bond!=NULL;bond=bond->next)
    {
      n=Get_Atom_Number(bond->atom);
      fprintf(MF,"%4d",n);
    }
    fprintf(MF,"\n");
  }
}

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

int Save_Steric_Data(FILE *MF, Ster *ster)
{
  int i=0,n=0;
  Conf *conf=NULL;
  if(ster==NULL) return(0);
  fprintf(MF,"STERIC  %c %10.6f%10.6f%10.6f%10.6f%10.6f%10.6f%10.6f  %s\n",ster->type
            ,ster->tot_val,ster->err_val
            ,ster->tot_con,ster->err_con
            ,ster->max_val,ster->pr_area
            ,ster->peak_R,ster->name);
  if(ster->val!=NULL)
  {
    fprintf(MF,"PROFIL  %10.6f%10.6f%6d  %s\n"
              ,ster->min,ster->max,ster->size,ster->name);
    for(i=0;isize;)
    {
      for(n=0;n<10;n++,i++)
      {
        fprintf(MF,"%10.6f",ster->val[i]);
      }
      fprintf(MF,"\n");
    }
  }
  if(ster->conf!=NULL)
  {
    fprintf(MF,"CONFS   %10.6f%10.6f  %s\n"
              ,ster->tot_con,ster->err_con,ster->name);
    for(conf=First_Conformer(ster->conf);conf!=NULL;conf=conf->next)
    {
      fprintf(MF,"CONFOR  %10.6f%10.3f%10.6f%10.6f\n"
                ,conf->value,conf->temperature,conf->Ec,conf->mol);
    }
  }
  return(1);
}

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

void Save_Steric_Molecule(Mol *M, char *fname)
{
  FILE *MF=NULL;
  Ster *ster=NULL;
  char line[LINELEN];
  if(!strlen(fname)) strcpy(fname,M->Fname);
  New_Extension(fname,".stc");
  sprintf(line,"Saving molecule to file %s",fname);
  Out_Message(line,O_NEWLN);
  if((MF=fopen(fname,"wt"))==NULL)
  {
    Error_Message(E_OPNOUT,"Save Steric Molecule");
  }
  else
  {
    fprintf(MF,"STERIC %s\n",M->name);
    Save_Steric_Atoms(MF,M);
    for(ster=First_Steric(M->ster);ster!=NULL;ster=ster->next)
      Save_Steric_Data(MF,ster);
    fclose(MF);
  }
}

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

void Read_Steric_Atom(char *line, Mol *M)
{
  int n=0;
  char name[LINELEN],type[LINELEN];
  Vector a={0.0,0.0,0.0};

  if(sscanf(line,"%s%s%lf%lf%lf%s",name,type,&a.x,&a.y,&a.z,name)==6)
  {
    name[ATM_LEN]=0;
    type[ATM_LEN]=0;
    n++;
    M->atoms=New_Atom(M->atoms);
    Initialize_Atom(M->atoms);
    if(!M->main_atom) M->main_atom=M->atoms;
    strcpy(M->atoms->name,name);
    strcpy(M->atoms->type,type);
    M->atoms->v=VequalV(a);
    M->atoms->fv=VequalV(a);
  }
}

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

void Read_Steric_Bond(char *line, Mol *M)
{
  int A,b[12],i=0,num=0;             /* maximum coordination number of 12 */
  Atms *atomA=NULL, *atomB=NULL;

  if((num=sscanf(line,"%d%d%d%d%d%d%d%d%d%d%d%d%d",&A
		     ,&b[0],&b[1],&b[2],&b[3]
		     ,&b[4],&b[5],&b[6],&b[7]
		     ,&b[8],&b[9],&b[10],&b[11]))>1)
  {
    atomA=Goto_Atom(M->atoms,A);
    if((atomA)&&(atomA==M->origin)) M->main_atom=Goto_Atom(M->atoms,b[0]);
    for(i=0;iatoms,b[i]);
      Add_New_Bond(atomA,atomB,SINGLE_B);
      if((atomB)&&(atomB==M->origin)) M->main_atom=atomA;
    }
  }
}

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

Mol *Load_Steric_Molecule(FILE *MF, Mol *M)
{
  char line[LINELEN],name[LINELEN];
  int num,n;
  Atms *at=NULL;

/* test that file is really OK                                            */

  if((num=count_lines(MF,"ATOM"))==0)
  {
    Error_Message(E_NOATOM,"Load Steric Molecule");
	return(NULL);
  }
  num+=count_lines(MF,"ORIGIN");
  if(count_lines(MF,"BOND")==0) Error_Message(E_NOBOND,"Load Steric Molecule");
  if((M=New_Molecule(M))==NULL) return(NULL);
  Initialize_Molecule(M);
  rewind(MF);
  if(get_next_line(MF,"STERIC",line)==0)
  {
	Error_Message(E_NOTITL,"Load Steric Molecule");
	n=Get_Molecule_Number(M);
	sprintf(M->name,"M%d",n);
  }
  else sscanf(line,"%s%s",name,M->name);
  M->num_atoms=num;
  sprintf(line,"Loading steric data %s ...",M->name);
  Out_Message(line,O_NEWLN);

  while((fgets(line,LINELEN,MF)!=NULL)&&(!feof(MF)))
  {
    if(sscanf(line,"%s",name)<1) continue;
    if(strcmp(name,"ATOM")==0) Read_Steric_Atom(line+6,M);
    if(strcmp(name,"ORIGIN")==0)
    {
      Read_Steric_Atom(line+6,M);
      M->origin=M->atoms;
    }
    if(strcmp(name,"BOND")==0) Read_Steric_Bond(line+6,M);
  }

  if(M->atoms==NULL) n=0;
  else
  {
    for(at=First_Atom(M->atoms),n=0;at!=NULL;at=at->next,n++);
  }
  if(n!=num)
  {
	Error_Message(E_BDATOM,"Load Steric Molecule");
	M->num_atoms=n;
  }

/* read in bonding table                                                  */

  return(M);

}

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

Mol *Load_Alchemy_Molecule(FILE *MF, Mol *M)
{
  char line[LINELEN],name[LINELEN];
  int Natoms,Nbonds,n,i;
  char arb1[20],arb2[20],arb3[20];
  Vector a={0.0,0.0,0.0};
  int A,B;
  Atms *atomA=NULL, *atomB=NULL;

/* test that file is really OK                                            */

  if(fgets(line,98,MF)==NULL) return(NULL);
  if((n=sscanf(line,"%d%s%d%s%d%s%s"
     ,&Natoms,arb1,&Nbonds,arb2,&i,arb3,name))<3)
  {
    Error_Message(E_BDFILE,"Load Alchemy Molecule");
    return(NULL);
  }
  if((M=New_Molecule(M))==NULL) return(NULL);
  Initialize_Molecule(M);
  if(n<7)
  {
	Error_Message(E_NOTITL,"Load Alchemy Molecule");
	n=Get_Molecule_Number(M);
	sprintf(M->name,"M%d",n);
  }
  else strcpy(M->name,name);
  M->num_atoms=Natoms;
  sprintf(line,"Loading alchemy data %s ...",M->name);
  Out_Message(line,O_NEWLN);

/* read in atomic information                                             */

  n=0; i=0;
  while((nnum_atoms)&&(inum_atoms)&&(fgets(line,LINELEN,MF)!=NULL))
  {
    n++;
    M->atoms=New_Atom(M->atoms);
    Initialize_Atom(M->atoms);
    if(!M->main_atom) M->main_atom=M->atoms;
    sscanf(line,"%d%s%lf%lf%lf",&i,name,&a.x,&a.y,&a.z);
    name[ATM_LEN]=0;
    if (stricmp(name,"DUO")==0) M->origin=M->atoms;
    strcpy(M->atoms->name,name);
    strcpy(M->atoms->type,name);
    M->atoms->v=VequalV(a);
    M->atoms->fv=VequalV(a);
  }
  if(n!=i)
  {
	Error_Message(E_BDATOM,"Load Alchemy Molecule");
	M->num_atoms=Natoms;
  }

/* read in bonding table                                                  */

  n=0; i=0;
  while((natoms,A);
      atomB=Goto_Atom(M->atoms,B);
      if((atomA)&&(atomA==M->origin)) M->main_atom=atomB;
      if((atomB)&&(atomB==M->origin)) M->main_atom=atomA;
      if(strcmp(name,"SINGLE")==0) Add_New_Bond(atomA,atomB,SINGLE_B);
      else if(strcmp(name,"DOUBLE")==0) Add_New_Bond(atomA,atomB,DOUBLE_B);
      else if(strcmp(name,"TRIPLE")==0) Add_New_Bond(atomA,atomB,TRIPLE_B);
      else if(strcmp(name,"AROMATIC")==0) Add_New_Bond(atomA,atomB,AROMAT_B);
      else Add_New_Bond(atomA,atomB,UNKNOW_B);
    }
  }

  return(M);
}

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

Mol *Load_Biograph_Molecule(FILE *MF, Mol *M)
{
  char line[LINELEN],name[LINELEN];
  int num,n;
  char arb1[LINELEN];
  Vector a={0.0,0.0,0.0};

/* test that file is really OK                                            */

  if((num=count_lines(MF,"HETATM"))==0)
  {
    Error_Message(E_NOATOM,"Load Biograph Molecule");
    return(NULL);
  }
  if(count_lines(MF,"CONECT")==0) Error_Message(E_NOBOND,"Load Biograph Molecule");
  if((M=New_Molecule(M))==NULL) return(NULL);
  Initialize_Molecule(M);
  if(get_next_line(MF,"DESCRP",line)==0)
  {
    Error_Message(E_NOTITL,"Load Biograph Molecule");
    n=Get_Molecule_Number(M);
    sprintf(M->name,"M%d",n);
  }
  else sscanf(line,"%s%s",arb1,M->name);
  M->num_atoms=num;
  sprintf(line,"Loading biograph data %s ...",M->name);
  Out_Message(line,O_NEWLN);

/* read in atomic information                                             */

  if(get_next_line(MF,"HETATM",line)==0)
  {
    Error_Message(E_NOATOM,"Load Biograph Molecule");
    return(Close_Current_Molecule(M));
  }
  n=0;
  do
  {
    line[6]=0;
    if(strcmp(line,"HETATM")!=0) break;
    if(sscanf(line+12,"%s",arb1)<1) break;
    arb1[ATM_LEN]=0;
    if(sscanf(line+BGF_SPACE+1,"%lf%lf%lf%s",&a.x,&a.y,&a.z,name)<4) continue;
    name[ATM_LEN]=0;
    n++;
    M->atoms=New_Atom(M->atoms);
    Initialize_Atom(M->atoms);
    if(!M->main_atom) M->main_atom=M->atoms;
    if (stricmp(name,"DUO")==0) M->origin=M->atoms;
    strcpy(M->atoms->name,arb1);
    strcpy(M->atoms->type,name);
    M->atoms->v=VequalV(a);
    M->atoms->fv=VequalV(a);
  } while(fgets(line,LINELEN,MF)!=NULL);
  if(n!=num)
  {
    Error_Message(E_BDATOM,"Load Biograph Molecule");
    M->num_atoms=n;
  }

/* read in bonding table                                                  */

  if(get_next_line(MF,"CONECT",line)==0)
  {
    Error_Message(E_NOBOND,"Load Biograph Molecule");
  }
  do
  {
    line[6]=0;
    if(strcmp(line,"CONECT")!=0) continue;
    Read_Steric_Bond(line+7,M);
  } while(fgets(line,98,MF)!=NULL);

  return(M);
}

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

Mol *Load_Biosym_Molecule(FILE *MF, Mol *M, char *file, Set *set)
{
  char line[LINELEN],name[LINELEN];
  int num,n,p;
  char arb1[LINELEN], arb2[LINELEN], arb3[LINELEN];
  Vector a={0.0,0.0,0.0};

  if((M=New_Molecule(M))==NULL) return(NULL);
  Initialize_Molecule(M);
  if(strrchr(file,'.')!=NULL) strrchr(file,'.')[0]=0;
  if(strlen(file)<1)
  {
    Error_Message(E_NOTITL,"Load Biosym Molecule");
    n=Get_Molecule_Number(M);
    sprintf(M->name,"M%d",n);
  }
  else strcpy(M->name,file);
  sprintf(line,"Loading biosymm data %s ...",M->name);
  Out_Message(line,O_NEWLN);

/* read in atomic information                                             */

  if(get_next_line(MF,"!DATE",line)==0)
  {
    Error_Message(E_NOATOM,"Load Biosym Molecule");
    return(Close_Current_Molecule(M));
  }
  n=0;
  while((fgets(line,LINELEN,MF)!=NULL)&&(!feof(MF)))
  {
    if(line[0]=='!') break;
    if(sscanf(line,"%s",name)<1) continue;
    strupr(name);
    if(strcmp(name,"END")==0) break;
    if((num=sscanf(line,"%s%lf%lf%lf%s%d%s%s",name,&a.x,&a.y,&a.z,arb3,&p,arb2,arb1))<6) continue;
    if(num<8)
    {
      strcpy(arb1,name);
      Screen_Atom_Name(arb1);
    }
    arb1[ATM_LEN]=0;
    name[ATM_LEN]=0;
    n++;
    M->atoms=New_Atom(M->atoms);
    Initialize_Atom(M->atoms);
    if(!M->main_atom) M->main_atom=M->atoms;
    if (stricmp(name,"DUO")==0) M->origin=M->atoms;
    strcpy(M->atoms->name,name);
    strcpy(M->atoms->type,arb1);
    M->atoms->v=VequalV(a);
    M->atoms->fv=VequalV(a);
  }
  M->num_atoms=n;

  Find_Atomic_Radii(M,set);
  Find_All_Bonds(M,set,0,0);

  return(M);
}

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

Mol *Load_PDB_Molecule(FILE *MF, Mol *M, char *file)
{
  char line[LINELEN],name[LINELEN];
  int num,n,p,q;
  char arb1[LINELEN];
  Vector a={0.0,0.0,0.0};
  Symm O,S;
  Atms *A=NULL;

  O.M.x=Vequal(1.0,0.0,0.0);
  O.M.y=Vequal(0.0,1.0,0.0);
  O.M.z=Vequal(0.0,0.0,1.0);
  O.t  =Vequal(0.0,0.0,0.0);

  S.M.x=Vequal(1.0,0.0,0.0);
  S.M.y=Vequal(0.0,1.0,0.0);
  S.M.z=Vequal(0.0,0.0,1.0);
  S.t  =Vequal(0.0,0.0,0.0);

/* test that file is really OK                                            */

  if((num=count_lines(MF,"ATOM"))==0)
  {
    Error_Message(E_NOATOM,"Load PDB Molecule");
    return(NULL);
  }
  if(count_lines(MF,"CONECT")==0) Error_Message(E_NOBOND,"Load PDB Molecule");
  if((M=New_Molecule(M))==NULL) return(NULL);
  Initialize_Molecule(M);
  if(strrchr(file,'.')!=NULL) strrchr(file,'.')[0]=0;
  if(strlen(file)<1)
  {
    Error_Message(E_NOTITL,"Load PDB Molecule");
    n=Get_Molecule_Number(M);
    sprintf(M->name,"M%d",n);
  }
  else strcpy(M->name,file);
  M->num_atoms=num;
  sprintf(line,"Loading PDB data %s ...",M->name);
  Out_Message(line,O_NEWLN);

/* read in atomic information                                             */

  n=0;
  while((fgets(line,LINELEN,MF)!=NULL)&&(!feof(MF)))
  {
    if(sscanf(line,"%s",name)<1) continue;
    strupr(name);
    if(strcmp(name,"ATOM")==0)
    {
      if(sscanf(line+4,"%d%s%d%lf%lf%lf",&p,arb1,&q,&a.x,&a.y,&a.z)<6) continue;
      sprintf(name,"%s%d",arb1,p);
      arb1[ATM_LEN]=0;
      name[ATM_LEN]=0;
      n++;
      M->atoms=New_Atom(M->atoms);
      Initialize_Atom(M->atoms);
      if(!M->main_atom) M->main_atom=M->atoms;
      if (stricmp(name,"DUO")==0) M->origin=M->atoms;
      strcpy(M->atoms->name,name);
      strcpy(M->atoms->type,arb1);
      M->atoms->v=VequalV(a);
      M->atoms->fv=VequalV(a);
    }
    else if(strcmp(name,"CONECT")==0) Read_Steric_Bond(line+6,M);
    else if(strcmp(name,"ORIGX1")==0) sscanf(line+6,"%lf%lf%lf%lf",&O.M.x.x,&O.M.x.y,&O.M.x.z,&O.t.x);
    else if(strcmp(name,"ORIGX2")==0) sscanf(line+6,"%lf%lf%lf%lf",&O.M.y.x,&O.M.y.y,&O.M.y.z,&O.t.y);
    else if(strcmp(name,"ORIGX3")==0) sscanf(line+6,"%lf%lf%lf%lf",&O.M.z.x,&O.M.z.y,&O.M.z.z,&O.t.z);
    else if(strcmp(name,"SCALE1")==0) sscanf(line+6,"%lf%lf%lf%lf",&S.M.x.x,&S.M.x.y,&S.M.x.z,&S.t.x);
    else if(strcmp(name,"SCALE2")==0) sscanf(line+6,"%lf%lf%lf%lf",&S.M.y.x,&S.M.y.y,&S.M.y.z,&S.t.y);
    else if(strcmp(name,"SCALE3")==0) sscanf(line+6,"%lf%lf%lf%lf",&S.M.z.x,&S.M.z.y,&S.M.z.z,&S.t.z);
    else if(strcmp(name,"END")==0) break;
  }
  
  if(n!=num)
  {
    Error_Message(E_BDATOM,"Load PDB Molecule");
    M->num_atoms=n;
  }

  if(!Identity_Operator(&O))
  {
    Perform_Symmetry(M,&O);
    for(A=First_Atom(M->atoms);A!=NULL;A=A->next) A->v=VequalV(A->fv);
  }

  return(M);
}

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

Mol *Load_Schakal_Molecule(FILE *MF, Mol *M)
{
  char line[LINELEN],name[LINELEN];
  int num,n;
  char arb1[LINELEN];

/* test that file is really OK                                            */

  if((num=count_lines(MF,"AT"))==0)
  {
    Error_Message(E_NOATOM,"Load Schakal Molecule");
    return(NULL);
  }
  if((M=New_Molecule(M))==NULL) return(NULL);
  Initialize_Molecule(M);
  M->mode|=FRAC_BIT;              /* default to expect fractional coords. */ 

  if(get_next_line(MF,"TITL",line)==0)
  {
    Error_Message(E_NOTITL,"Load Schakal Molecule");
    n=Get_Molecule_Number(M);
    sprintf(M->name,"M%d",n);
  }
  else sscanf(line,"%s%s",arb1,M->name);
  M->num_atoms=num;
  sprintf(line,"Loading Schakal data %s ...",M->name);
  Out_Message(line,O_NEWLN);

/* read in unit cell information, if any                                  */

  if(get_next_line(MF,"CE",line)==0)
  {
    Error_Message(E_NOCELL,"Load Schakal Molecule");
    M->mode&=(FULL_BIT^FRAC_BIT);
  }
  else Get_Cell_Parameters(M,strchr(line,' '),F_SCHAKAL);

/* read in all other information (atomic, symmetry)                       */

  n=0;
  while((fgets(line,LINELEN,MF)!=NULL)&&(!feof(MF)))
  {
    if(sscanf(line,"%s",name)<1) continue;
    name[2]=0;
    strupr(name);
    if(strcmp(name,"AT")==0) n+=Get_Atom_Coordinates(M,strchr(line,' '),F_SCHAKAL);
    else if(strcmp(name,"SY")==0) Get_Symmetry_Operator(M,strchr(line,' '),F_SCHAKAL,S_NORM);
    else if(strcmp(name,"DU")==0) Get_Symmetry_Operator(M,strchr(line,' '),F_SCHAKAL,S_CENT);
  }

  if(n!=num)
  {
    Error_Message(E_BDATOM,"Load Schakal Molecule");
    M->num_atoms=n;
  }
  if(M->mode&CSSG_BIT) M->geneq*=2;

  return(M);
}

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

Mol *Load_Shelx_Molecule(FILE *MF, Mol *M)
{
  char line[LINELEN],name[LINELEN];
  int n;
  char arb1[LINELEN];

  if((M=New_Molecule(M))==NULL) return(NULL);
  Initialize_Molecule(M);
  M->mode|=FRAC_BIT;              /* default to expect fractional coords. */ 

  if(get_next_line(MF,"TITL",line)==0)
  {
    Error_Message(E_NOTITL,"Load Shelx Molecule");
    n=Get_Molecule_Number(M);
    sprintf(M->name,"M%d",n);
  }
  else sscanf(line,"%s%s",arb1,M->name);
  sprintf(line,"Loading shelx data %s ...",M->name);
  Out_Message(line,O_NEWLN);

/* read in all other information (atomic, symmetry)                       */

  n=0;
  while((n>=0)&&(fgets(line,LINELEN,MF)!=NULL)&&(!feof(MF)))
  {
    if((line[0]==' ')||(sscanf(line,"%s",name)!=1)) continue;
    strupr(name);
    switch(Shelx_Card_Type(name))
    {
      case END : n*=-1; break;
      case CEL : Get_Cell_Parameters(M,strchr(line,' '),F_SHELXL); break;
      case LAT : Get_Lattice(M,strchr(line,' '),F_SHELXL); break;
      case SYM : Get_Symmetry_Operator(M,strchr(line,' '),F_SHELXL,S_NORM); break;
      case ATM : n+=Get_Atom_Coordinates(M,line,F_SHELXL); break;
      default  : break;
    }
  }
  if(n<0) n*=-1;
  M->num_atoms=n;
  if(M->mode&CSSG_BIT) M->geneq*=2;

  return(M);
}

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

Mol *Load_Gstcor_Molecule(FILE *MF, Mol *M, Set *set)
{
  char line[LINELEN],name[LINELEN];
  int n;

  while(fgets(line,LINELEN,MF),!feof(MF))
  {
    if(strncmp(line+8,"**FRAG**",8)==0) break;
  }

  while(strncmp(line+8,"**FRAG**",8)==0)
  {
    if((M=New_Molecule(M))==NULL) return(NULL);
    Initialize_Molecule(M);
    M->mode|=FRAC_BIT;            /* default to expect fractional coords. */ 
    line[8]=0;
    n=strlen(line);
    while(n--,line[n]==' ') line[n]=0;
    strcpy(M->name,line);
    sprintf(line,"Loading gstat coordinate data %s ...",M->name);
    Out_Message(line,O_NEWLN);

/* read in unit cell information, if any                                  */

    n=0;
    while((fgets(line,LINELEN,MF)!=NULL)&&(!feof(MF)))
    {
      if(strncmp(line+8,"**FRAG**",8)==0) break;
      if(sscanf(line,"%s",name)<1) continue;
      strupr(name);
      if(strncmp(name,"SYMM",4)==0) Get_Symmetry_Operator(M,strchr(line,' '),F_GSTCOR,S_NORM);
      else if(strncmp(name,"CELL",4)==0) Get_Cell_Parameters(M,strchr(line,' '),F_GSTCOR);
      else n+=Get_Atom_Coordinates(M,line,F_GSTCOR);
    }
    M->num_atoms=n;
    if(M->mode&CSSG_BIT) M->geneq*=2;

/* create cartesian coordinate vectors                                    */

    Create_Unit_Cell_Atoms(M);
    Convert_to_Cartesian(M);
    Setup_Atomic_Parameters(M,set);
    Find_All_Bonds(M,set,0,0);
    Find_All_Bonded_Groups(M,0);

  }

  return(M);
}

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

Mol *Load_Xtal_Molecule(FILE *MF, Mol *M)
{
  char line[LINELEN],name[LINELEN];
  int num,n;
  char arb1[LINELEN];

/* test that file is really OK                                            */

  if((num=count_lines(MF,"ATOM"))==0)
  {
    Error_Message(E_NOATOM,"Load Xtal Molecule");
    return(NULL);
  }
  if((M=New_Molecule(M))==NULL) return(NULL);
  Initialize_Molecule(M);
  M->mode|=FRAC_BIT;              /* default to expect fractional coords. */ 

  if(get_next_line(MF,"TITL",line)==0)
  {
    Error_Message(E_NOTITL,"Load Xtal Molecule");
    n=Get_Molecule_Number(M);
    sprintf(M->name,"M%d",n);
  }
  else sscanf(line,"%s%s",arb1,M->name);
  M->num_atoms=num;
  sprintf(line,"Loading xtal data %s ...",M->name);
  Out_Message(line,O_NEWLN);

/* read in all other information (atomic, symmetry)                       */

  n=0;
  while((fgets(line,LINELEN,MF)!=NULL)&&(!feof(MF)))
  {
    if(sscanf(line,"%s",name)<1) continue;
    strupr(name);
    if(strcmp(name,"ATOM")==0) n+=Get_Atom_Coordinates(M,strchr(line,' '),F_XTAL);
    else if(strcmp(name,"SYMTRY")==0) Get_Symmetry_Operator(M,strchr(line,' '),F_XTAL,S_NORM);
    else if(strcmp(name,"LATICE")==0) Get_Lattice(M,strchr(line,' '),F_XTAL);
    else if(strcmp(name,"CELL")==0) Get_Cell_Parameters(M,strchr(line,' '),F_XTAL);
  }

  if(n!=num)
  {
    Error_Message(E_BDATOM,"Load Xtal Molecule");
    M->num_atoms=n;
  }
  if(M->mode&CSSG_BIT) M->geneq*=2;

  return(M);
}

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

Mol *Load_Acryst_Molecule(FILE *MF, Mol *M)
{
  char line[LINELEN];
  int num,n;

/* test that file is really OK                                            */

  if(fgets(line,LINELEN,MF)==NULL) return(NULL);
  if((M=New_Molecule(M))==NULL) return(NULL);
  Initialize_Molecule(M);
  M->mode|=FRAC_BIT;              /* default to expect fractional coords. */ 

  if(strlen(line)<1)
  {
    Error_Message(E_NOTITL,"Load Acryst Molecule");
    n=Get_Molecule_Number(M);
    sprintf(M->name,"M%d",n);
  }
  else sscanf(line,"%s",M->name);
  sprintf(line,"Loading alchemy crystal data %s ...",M->name);
  Out_Message(line,O_NEWLN);

  if(fgets(line,LINELEN,MF)==NULL) return(Close_Current_Molecule(M));

  if((sscanf(line,"%d",&num)!=1)||(num<1))
  {
    Error_Message(E_NOATOM,"Load Acryst Molecule");
    return(Close_Current_Molecule(M));
  }
  M->num_atoms=num;

  if(fgets(line,LINELEN,MF)==NULL)
  {
    Error_Message(E_NOCELL,"Load Acryst Molecule");
    return(Close_Current_Molecule(M));
  }
  if(!Get_Cell_Parameters(M,line,F_ACRYST)) Error_Message(E_BADCEL,"Load Acryst Molecule");

  for(fgets(line,LINELEN,MF),n=0;nnum_atoms=n;
  }
  if(M->mode&CSSG_BIT) M->geneq*=2;

  return(M);
}

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

void New_Extension(char *file, char *newext)
{
  char *ext;
  ext=strrchr(file,'.');
  if((ext!=NULL)&&(ext[-1]=='.')) ext=NULL;
  if(ext!=NULL) ext[0]=0;
  if(newext!=NULL) strcat(file,newext);
}

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

unsigned Get_File_Mode(FILE *in)
{
  char line[LINELEN];
  char line2[LINELEN];
  char line3[LINELEN];
  char f[6][20];
  int num=0;

  rewind(in);
  if(fgets(line,LINELEN,in)==NULL)
    if(feof(in)) return(0);
  if(fgets(line2,LINELEN,in)==NULL)
    if(feof(in)) return(0);
  if(fgets(line3,LINELEN,in)==NULL)
    if(feof(in)) return(0);
  rewind(in);
  strupr(line);
  strupr(line2);
  strupr(line3);
  sscanf(line,"%s%s%s%s%s%s",f[0],f[1],f[2],f[3],f[4],f[5]);

  if(strncmp(line,"#STERIC",7)==0) return(F_STERIN);
  if(strncmp(line,"#STERPAR",8)==0) return(F_STERPAR);
  if(strncmp(line,"STERIC",6)==0) return(F_STERIC);
  if(strncmp(line,"BIOGRF",6)==0) return(F_BIOGRF);
  if(strncmp(line,"!BIOSYM",7)==0) return(F_BIOSYM);
  if((strncmp(line ,"ORIGX1",6)==0)
   &&(strncmp(line2,"ORIGX2",6)==0)
   &&(strncmp(line3,"ORIGX3",6)==0)) return(F_BIOPDB);
  if((strncmp(line,"TITL",4)==0)&&(strncmp(line2,"CELL",4)==0))
  {
    if(strncmp(line3,"ZERR",4)==0) return(F_SHELXL);
    if(strncmp(line3,"LATT",4)==0) return(F_SHELXL);
    return(F_SCHAKAL);
  }
  if(strncmp(line+8,"**FRAG**",8)==0) return(F_GSTCOR);
  if((strcmp(f[1],"ATOMS,")==0)
   &&(strcmp(f[3],"BONDS,")==0)
   &&(strncmp(f[5],"CHARGES,",7)==0)) return(F_ALCHEM);
  if(strncmp(line,"HEADER:",7)==0) return(F_CONFOR);
  if(strncmp(line2,"COMPID",6)==0) return(F_XTAL);
  if((strncmp(line2,"   ",3)==0)&&(sscanf(line2,"%d",&num)==1)&&(num>0)) return(F_ACRYST);
  return(0);
}

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

FILE *Open_Input_File(char *file, unsigned *mode)
{
  FILE *in=NULL;

  if((in=fopen(file,"rt"))==NULL) {
   if(New_Extension(file,".inp"),(in=fopen(file,"rt"))==NULL) {
    if(New_Extension(file,".stc"),(in=fopen(file,"rt"))==NULL) {
     if(New_Extension(file,".bgf"),(in=fopen(file,"rt"))==NULL) {
      if(New_Extension(file,".mol"),(in=fopen(file,"rt"))==NULL) {
       if(New_Extension(file,".trj"),(in=fopen(file,"rt"))==NULL) {
        if(New_Extension(file,".dat"),(in=fopen(file,"rt"))==NULL) {
         if(New_Extension(file,".cor"),(in=fopen(file,"rt"))==NULL) {
          if(New_Extension(file,".cry"),(in=fopen(file,"rt"))==NULL) {
           if(New_Extension(file,".shl"),(in=fopen(file,"rt"))==NULL) {
            if(New_Extension(file,".res"),(in=fopen(file,"rt"))==NULL) {
             if(New_Extension(file,".car"),(in=fopen(file,"rt"))==NULL) {
              if(New_Extension(file,".pdb"),(in=fopen(file,"rt"))==NULL) {
               if(New_Extension(file,""),(in=fopen(file,"rt"))==NULL) {
                Error_Message(E_BDFILE,"Open Input File");
               }
              }
             }
            }
           }
          }
         }
        }
       }
      }
     }
    }
   }
  }
  Out_Message("Opening file ",O_BLANK);
  Out_Message(file,O_NEWLN);
  *mode=Get_File_Mode(in);
  return(in);
}

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

Mol *Load_New_Molecule(char *file, Mol *current, Set *set)
{
  FILE *in;
  unsigned mode=F_STERIC;
  Mol *new=NULL;

  if((in=Open_Input_File(file,&mode))==NULL) return(NULL);
  switch(mode)
  {
    case F_STERPAR: Load_Parameters(set,in,file); return(current); break;
    case F_STERIN : if(set->input!=stdin) fclose(set->input);
                    set->input=in;
                    return(new); break;
    case F_STERIC : new=Load_Steric_Molecule(in,current);
		            break;
    case F_BIOGRF : new=Load_Biograph_Molecule(in,current);
		    break;
    case F_ALCHEM : new=Load_Alchemy_Molecule(in,current);
		    break;
    case F_BIOSYM : new=Load_Biosym_Molecule(in,current,file,set);
		    break;
    case F_BIOPDB : new=Load_PDB_Molecule(in,current,file);
		    break;
    case F_SCHAKAL: new=Load_Schakal_Molecule(in,current);
		    break;
    case F_SHELXL : new=Load_Shelx_Molecule(in,current);
		    break;
    case F_GSTCOR : new=Load_Gstcor_Molecule(in,current,set);
		    break;
    case F_ACRYST : new=Load_Acryst_Molecule(in,current);
		    break;
    case F_XTAL   : new=Load_Xtal_Molecule(in,current);
		    break;
    case F_CONFOR : Error_Message(E_OPCONF,"Load New Molecule"); break;
    default       : Error_Message(E_FTYPE,"Load New Molecule"); break;
  }
  fclose(in);
  if(new==NULL) Error_Message(E_LOAD,"Load New Molecule");
  else
  {
    if((mode==F_SCHAKAL)||(mode==F_SHELXL)||(mode==F_ACRYST)||(mode==F_XTAL))
    {
      Create_Unit_Cell_Atoms(new);
      Convert_to_Cartesian(new);
      Setup_Atomic_Parameters(new,set);
      Find_All_Bonds(new,set,0,0);
      Find_All_Bonded_Groups(new,0);
    }
    else if(mode!=F_GSTCOR) Setup_Atomic_Parameters(new,set);
    New_Extension(file,"");
    strcpy(new->Fname,file);
  }
  return(new);
}

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