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

/**************************************************************************/
/******************  Atom 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 "crystal.h"       /* functions for all file manipulations        */
#include "stermem.h"       /* dynamic memory management module            */
#include "stercalc.h"      /* main calculation routines                   */
#include "sterfile.h"      /* functions for all file manipulations        */
#include "stertext.h"      /* all text functions for text mode            */
#include "stergrp.h"       /* group manipulation functions                */
#include "steraid.h"       /* additional functions needed                 */

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

int Initialize_Symmetry(Symm *symm, char *S, Matrix M, Vector t, unsigned type)
{
  strncpy(symm->S,S,SYM_LEN);
  symm->M.x=VequalV(M.x);
  symm->M.y=VequalV(M.y);
  symm->M.z=VequalV(M.z);
  symm->t=VequalV(t);
  symm->mode=0;
  switch(type)
  {
    case S_NORM : symm->mode=0; break;
    case S_CENT : symm->mode|=CENT_BIT; break;
    default     : break;
  }
  return(1);
}

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

int Copy_Atom(Atms *A, Atms *B)
{
  strcpy(A->name,B->name);
  strcpy(A->type,B->type);
  A->v=VequalV(B->v);
  A->fv=VequalV(B->fv);
  A->distance=B->distance;
  A->theta=B->theta;
  A->phi=B->phi;
  A->radius=B->radius;
  A->SVangle=B->SVangle;
  A->stat=B->stat;
  return(1);
/* next and prev pointers have already been assigned, do not fiddle !!!   */
}

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

int Identity_Operator(Symm *symm)
{
  if(Vcmp(symm->M.x,Vequal(1.0,0.0,0.0))) return(0);
  if(Vcmp(symm->M.y,Vequal(0.0,1.0,0.0))) return(0);
  if(Vcmp(symm->M.z,Vequal(0.0,0.0,1.0))) return(0);
  if(Vcmp(symm->t,  Vequal(0.0,0.0,0.0))) return(0);
  return(1);
}

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

int Convert_to_Cartesian(Mol *M)
{
  double cral, crbe, crga, srga;
  double r,p,pa,qa;
  double cab,cac,cbb,cbc,ccc;
  Atms *A=NULL;
  
  if(M==NULL)
  {
    Error_Message(E_NOMOL,"Convert_to_Cartesian");
    return(0);
  }
  cral=cos(Radians(M->al));
  crbe=cos(Radians(M->be));
  crga=cos(Radians(M->ga));
  srga=sin(Radians(M->ga));
  pa = crbe * crga;
  r = (cral - pa) / srga;
  p = crbe;
  qa = sqrt((double)1. - p * p - r * r);
  cab = M->b * crga;
  cac = M->c * crbe;
  cbb = M->b * srga;
  cbc = M->c * r;
  ccc = M->c * qa;
  for(A=First_Atom(M->atoms);A!=NULL;A=A->next)
  {
    A->v.x = A->fv.x * M->a + A->fv.y * cab + A->fv.z * cac;
    A->v.y = A->fv.y * cbb  + A->fv.z * cbc;
    A->v.z = A->fv.z * ccc;
  }
  return(1);
}

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

int Create_Unit_Cell_Atoms(Mol *M)
{
  Atms *atom=NULL;
  Atms *new=NULL;
  int i;

  if(M==NULL) return(0);
  atom=Last_Atom(M->atoms);
  for(i=0;i<8;i++)
  {
    if((new=New_Atom(atom))==NULL) return(0);
    atom=new;
    M->num_atoms++;
    Initialize_Atom(atom);
    sprintf(atom->name,"A%d",i);
    sprintf(atom->type,"A");
    atom->group=FULL_BIT;
    atom->stat|=GRP_BIT;
    switch(i)
    {
      case 0  : atom->fv=Vequal(0.0,0.0,0.0); break;
      case 1  : atom->fv=Vequal(1.0,0.0,0.0); break;
      case 2  : atom->fv=Vequal(0.0,1.0,0.0); break;
      case 3  : atom->fv=Vequal(0.0,0.0,1.0); break;
      case 4  : atom->fv=Vequal(0.0,1.0,1.0); break;
      case 5  : atom->fv=Vequal(1.0,0.0,1.0); break;
      case 6  : atom->fv=Vequal(1.0,1.0,0.0); break;
      case 7  : atom->fv=Vequal(1.0,1.0,1.0); break;
      default : break;
    }
  }
  return(1);
}

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

int Remove_Duplicate_Atoms(Mol *M, Set *set)
{
  Atms *A=NULL;
  Atms *B=NULL;
  Grps *group=NULL;
  int i,n;
  char line[LINELEN];
  
  if(M==NULL) return(0);
  if((A=First_Atom(M->atoms))==NULL) return(0);
  
  while(A!=NULL)
  {
    B=A->next;
    while(B!=NULL)
    {
      if((A!=B)
       &&((!(set->modev&SNAM_BIT))||(strcmp(A->name,B->name)==0))
       &&(strcmp(A->type,B->type)==0)
       &&(!Vcmp(A->fv,B->fv)))
      {
        sprintf(line,"Removing duplicate atom %s",A->name);
        Out_Message(line,O_NEWLN);
        B=Close_Atom(B);
      }
      if(B!=NULL) B=B->next;
    }
    A=A->next;
  }
  for(i=0,A=First_Atom(M->atoms);A!=NULL;A=A->next,i++);
  M->num_atoms=i;
  for(group=First_Group(M->groups);group!=NULL;group=group->next)
  {
    n=Get_Group_Number(group);
    for(i=0,A=First_Atom(M->atoms);A!=NULL;A=A->next)
    {
      if(A->group==n) i++;
    }
    group->num_atoms=i;
  }
  return(1);
}

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

/* void Add_New_Bond(Atms *atomA, Atms *atomB, short unsigned type) */
/* double Get_Atomic_Radius(int bonds, char *name, Parm *parm, unsigned mode) */

int Find_All_Bonds(Mol *M, Set *set, int gnum, unsigned mode)
{
  Atms *A=NULL;
  Atms *B=NULL;
  double rA,rB,rD;
  char line[LINELEN];
  int num=0;

  if(M==NULL) return(0);
  if(M->atoms==NULL) return(0);
  for(A=First_Atom(M->atoms);A!=NULL;A=A->next)
  {
    if((gnum!=0)&&(A->group!=gnum)) continue;
    rA=Get_Atomic_Radius(0,A->type,set->parm,CVL_BIT);
    for(B=A->next;B!=NULL;B=B->next)
    {
      if((gnum!=0)&&(B->group!=gnum)) continue;
      rB=Get_Atomic_Radius(0,B->type,set->parm,CVL_BIT);
      rD=Vmag(Vdif(A->v,B->v));
      if((((double)set->bond_minf*rA)+((double)set->bond_minf*rB)bond_maxf*rA)+((double)set->bond_maxf*rB)>rD))
      {
        Add_New_Bond(A,B,SINGLE_B);
        if(mode&VIS_BIT)
        {
          sprintf(line,"Atom %s (r=%4.2f) bonded to atom %s (r=%4.2f), with %8.4f seperation)"
                      ,A->name,rA,B->name,rB,rD);
          Out_Message(line,O_NEWLN);
        }
        num++;
      }
    }
  }
  if(num)
  {
    sprintf(line,"%d individual bonds were found in %s",num,M->name);
    Out_Message(line,O_NEWLN);
  }
  return(1);
}

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

int Perform_Symmetry(Mol *M, Symm *S)
{
  Atms *A=NULL;
  if(M==NULL) return(0);
  if(S==NULL) return(0);
  for(A=First_Atom(M->atoms);A!=NULL;A=A->next)
  {
    A->fv=Mtransform(S->M,A->fv);
    A->fv=Vsum(S->t,A->fv);
  }
  return(1);
}

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

Atms *Merge_Atom_Lists(Atms *A, Atms *B)
{
  if(B==NULL) return(First_Atom(A));
  if(A==NULL) return(First_Atom(B));
  if(((A=Last_Atom(A))!=NULL)&&((B=First_Atom(B))!=NULL))
  {
    B->prev=A;
    A->next=B;
  }
  return(First_Atom(A));
}

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

Atms *New_Atom_Symmetry(Mol *M, Symm *S, Atms *B, char *string)
{
  Atms *A=NULL;
  char line[LINELEN];

  if(M==NULL) return(B);
  if(S==NULL) return(B);
  if(M->atoms==NULL) return(B);
  sprintf(line,"%sAdding atoms with symmetry:",string); Out_Message(line,O_NEWLN);
  sprintf(line,"%5.1f%5.1f%5.1f%8.3f",S->M.x.x,S->M.x.y,S->M.x.z,S->t.x); Out_Message(line,O_NEWLN);
  sprintf(line,"%5.1f%5.1f%5.1f%8.3f",S->M.y.x,S->M.y.y,S->M.y.z,S->t.y); Out_Message(line,O_NEWLN);
  sprintf(line,"%5.1f%5.1f%5.1f%8.3f",S->M.z.x,S->M.z.y,S->M.z.z,S->t.z); Out_Message(line,O_NEWLN);
  for(A=First_Atom(M->atoms);A!=NULL;A=A->next)
  {
    if(strcmp(A->type,"A")==0) continue;
    if((B=New_Atom(B))==NULL) break;
    Initialize_Atom(B);
    Copy_Atom(B,A);
    B->fv=Mtransform(S->M,B->fv);
    B->fv=Vsum(S->t,B->fv);
  }
  return(B);
}

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

int Get_Symmetry_Matrix(char *S, Matrix *M, Vector *t)
{
  int i,len;
  M->x=Vequal(1.0,0.0,0.0);
  M->y=Vequal(0.0,1.0,0.0);
  M->z=Vequal(0.0,0.0,1.0);
  *t=Vequal(0.0,0.0,0.0);
  len=strlen(S);
  for(i=0;i'9'))) return(0);
  if(sscanf(S,"%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf"
            ,&M->x.x,&M->x.y,&M->x.z,&t->x
            ,&M->y.x,&M->y.y,&M->y.z,&t->y
            ,&M->z.x,&M->z.y,&M->z.z,&t->z)<12) return(0);
  return(1);
}

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

int do_M(Matrix *M, int row, int col, int val)
{
  if(row==2)
  {
    if(col==2) M->z.z=val;
    else if(col==1) M->z.y=val;
    else M->z.x=val;
  }
  else if(row==1)
  {
    if(col==2) M->y.z=val;
    else if(col==1) M->y.y=val;
    else M->y.x=val;
  }
  else
  {
    if(col==2) M->x.z=val;
    else if(col==1) M->x.y=val;
    else M->x.x=val;
  }
  return(1);
}

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

double Get_Tval(char *S, int *i)
{
  double value=0.0;
  double div=1.0;
  char str[LINELEN];
  int n;
  
  for(n=*i;(((S[n]=='+')&&(S[n+1]<=57)&&(S[n+1]>=46))
          ||((S[n]=='-')&&(S[n+1]<=57)&&(S[n+1]>=46))
          ||((S[n]<=57)&&(S[n]>=46)));n++);
  strncpy(str,S+*i,n-*i);
  str[n-*i]=0;
  sscanf(str,"%lf",&value);
  if(strchr(str,'/')!=NULL) sscanf(strchr(str,'/')+1,"%lf",&div);
  value/=div;
  *i=n-1;
  return(value);
}

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

int Get_Symmetry(char *S, Matrix *M, Vector *t)
{
  int i,len,p=0;
  double T[3]={0.0,0.0,0.0};
  len=strlen(S);
  strupr(S);
  M->x=Vequal(0.0,0.0,0.0);
  M->y=Vequal(0.0,0.0,0.0);
  M->z=Vequal(0.0,0.0,0.0);
  *t=Vequal(0.0,0.0,0.0);
  for(i=0;i=47))
     ||((S[i]=='-')&&(S[i+1]<=57)&&(S[i+1]>=47))
     ||((S[i]<=57)&&(S[i]>=47))) T[p]=Get_Tval(S,&i);
    else if((S[i]=='X')&&(S[i-1]!='-')) do_M(M,p,0,1);
    else if((S[i]=='-')&&(S[i+1]=='X')) do_M(M,p,0,-1);
    else if((S[i]=='Y')&&(S[i-1]!='-')) do_M(M,p,1,1);
    else if((S[i]=='-')&&(S[i+1]=='Y')) do_M(M,p,1,-1);
    else if((S[i]=='Z')&&(S[i-1]!='-')) do_M(M,p,2,1);
    else if((S[i]=='-')&&(S[i+1]=='Z')) do_M(M,p,2,-1);
  }
  *t=Vequal(T[0],T[1],T[2]);
  return(1);
}

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

int Transform_Data(Mol *M, char *SS, Set *set)
{
  Symm S;
  char line[LINELEN];

  if(M==NULL) return(0);
  if(M->atoms==NULL) return(0);

  if(strlen(SS)<2)
  {
    sprintf(SS,"Enter symmetry transform (x,y,z): ");
    if(Get_Input_Line(SS,set->input)==0) return(0);
    if(strlen(SS)<2) return(0);
  }
  strcpy(S.S,SS);
  
  if(!Get_Symmetry(S.S,&S.M,&S.t)) return(0);
  
  sprintf(line,"Symmetry Matrix and translation vector"); Out_Message(line,O_NEWLN);
  sprintf(line,"%5.1f%5.1f%5.1f%8.3f",S.M.x.x,S.M.x.y,S.M.x.z,S.t.x); Out_Message(line,O_NEWLN);
  sprintf(line,"%5.1f%5.1f%5.1f%8.3f",S.M.y.x,S.M.y.y,S.M.y.z,S.t.y); Out_Message(line,O_NEWLN);
  sprintf(line,"%5.1f%5.1f%5.1f%8.3f",S.M.z.x,S.M.z.y,S.M.z.z,S.t.z); Out_Message(line,O_NEWLN);
  
  Perform_Symmetry(M,&S);
  Convert_to_Cartesian(M);
  Calculate_Parameters(M,Vequal(0.0,0.0,0.0),set);
  return(1);
}

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

int Expand_Data(Mol *M, Set *set)
{
  int i;
  Symm *S=NULL;
  Atms *A=NULL;
  Symm C;
  char line[LINELEN];

  if(M==NULL) return(0);
  if(M->atoms==NULL) return(0);
  if(M->symmetry==NULL) return(0);

/* perform centrosymmetric symmetry operation                              */

  if(M->mode&CSSG_BIT)
  {
    C.M.x=Vequal(-1.0,0.0,0.0);
    C.M.y=Vequal(0.0,-1.0,0.0);
    C.M.z=Vequal(0.0,0.0,-1.0);
    C.t=Vequal(0.0,0.0,0.0);
    sprintf(line,"Centrosymmetric Operator - ");
    A=New_Atom_Symmetry(M,S,A,line);
    M->atoms=Merge_Atom_Lists(M->atoms,A);
    A=NULL;
  }

/* perform non centring symmetry operations                                */

  for(i=1,S=First_Symm(M->symmetry);S!=NULL;S=S->next,i++)
  {
    if(S->mode&CENT_BIT) continue;
    if(Identity_Operator(S)) continue;
    sprintf(line,"Symmetry operator %5d : ",i);
    A=New_Atom_Symmetry(M,S,A,line);
  }
  M->atoms=Merge_Atom_Lists(M->atoms,A);
  A=NULL;

/* perform centring symmetry operations (eg DU lines)                      */

  for(i=1,S=First_Symm(M->symmetry);S!=NULL;S=S->next,i++)
  {
    if(!(S->mode&CENT_BIT)) continue;
    if(Identity_Operator(S)) continue;
    sprintf(line,"Centring Condition %5d : ",i);
    A=New_Atom_Symmetry(M,S,A,line);
  }
  M->atoms=Merge_Atom_Lists(M->atoms,A);
  for(i=0,A=First_Atom(M->atoms);A!=NULL;A=A->next,i++);
  M->num_atoms=i;

/* remove now redundant symmetries                                         */

  M->symmetry=Close_All_Symmetries(M->symmetry);

  Remove_Duplicate_Atoms(M,set);

  Convert_to_Cartesian(M);
  Calculate_Parameters(M,Vequal(0.0,0.0,0.0),set);

  return(1);
}

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

int Box_Atoms(Mol *M, char *inline, Set *set)
{
  int i,n;
  char line[LINELEN];
  double a[3]={1.0,1.0,1.0};
  double m[3]={0.0,0.0,0.0};
  Symm S;
  Atms *A=NULL;
  Atms *B=NULL;

  if(M==NULL) return(0);
  if(M->atoms==NULL) return(0);

  n=sscanf(inline,"%lf%lf%lf%lf%lf%lf",&m[0],&m[1],&m[2],&a[0],&a[1],&a[2]);
  if(n<0) n=0;
  if(n==0) { m[0]=1.0; n=1; }
  if(n<3) m[2]=m[0];
  if(n<2) m[1]=m[0];
  if(n<4)
  {
    for(i=0;i<3;i++)
    {
      a[i]=m[i];
      if(m[i]<1.0) a[i]*=-1.0;
      else m[i]=0.0;
    }
  }
  
  sprintf(line,"Boxing from [%4.1f,%4.1f,%4.1f] to [%4.1f,%4.1f,%4.1f]"
              ,m[0],m[1],m[2],a[0],a[1],a[2]);
  Out_Message(line,O_NEWLN);

  S.S[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);

/* box into first unit cell                                                */

  Expand_Data(M,set);
  for(A=First_Atom(M->atoms);A!=NULL;A=A->next)
    A->fv=Vwrap(A->fv,Vequal(0.0,0.0,0.0),Vequal(1.0,1.0,1.0));

  if((a[0]>=1.0)||(a[1]>=1.0)||(a[2]>=1.0))
  {
    for(S.t.x=m[0];S.t.xatoms=Merge_Atom_Lists(M->atoms,B);
  for(i=0,A=First_Atom(M->atoms);A!=NULL;A=A->next,i++);
  M->num_atoms=i;

  if(set->modev&REDU_BIT) Remove_Duplicate_Atoms(M,set);

  Convert_to_Cartesian(M);
  Calculate_Parameters(M,Vequal(0.0,0.0,0.0),set);

  return(1);
}

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

char *Get_Type(char *name, char *type)
{
  int len;
  int i;
  strcpy(type,name);
  len=strlen(type);
  for(i=0;imode&=(FULL_BIT^FRAC_BIT);
  }
  else
  {
    M->a=a;
    M->b=b;
    M->c=c;
    M->al=al;
    M->be=be;
    M->ga=ga;
    cal=cos(Radians(al));
    cbe=cos(Radians(be));
    cga=cos(Radians(ga));
    M->cellvol=a*b*c*msqrt(1 - cal*cal - cbe*cbe - cga*cga + 2*cal*cbe*cga,"Get Cell Parameters : cell volume");
  }
  return(1);
}

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

int Get_Lattice(Mol *M, char *line, unsigned mode)
{
  int i=0;
  char centro[5]="c";
  char centri[5]="p";
  char xlat[9]="pirfabc";

  if(M==NULL) return(0);
  if(line==NULL) return(0);
  switch(mode)
  {
    case F_SHELXL : sscanf(line,"%d",&i);
                    if(i<0)
                    {
                      M->mode|=CSSG_BIT;
                      i*=-1;
                    }
                    if((i>0)&&(i<8)) M->centring=xlat[i-1];
                    break;
    case F_XTAL   : sscanf(line,"%s%s",centro,centri);
                    if((i=(int)(strchr(xlat,centri[0])-xlat))>=0)
                    {
                      M->centring=xlat[i];
                      if(centro[0]=='c') M->mode|=CSSG_BIT;
                    }
                    break;
    default       : break;
  }
  return(1);
}

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

int Get_Atom_Coordinates(Mol *M, char *line, unsigned mode)
{
  char name[LINELEN];
  char type[LINELEN];
  Vector a={0.0,0.0,0.0};
  int stype;

  if(M==NULL) return(0);
  if(line==NULL) return(0);

  switch(mode)
  {
    case F_SHELXL : if(sscanf(line,"%s%d%lf%lf%lf",name,&stype,&a.x,&a.y,&a.z)<5) return(0);
                    Get_Type(name,type);
                    break;
    case F_ACRYST : if(sscanf(line,"%s%s%lf%lf%lf",name,type,&a.x,&a.y,&a.z)<5) return(0);
                    break;
    default       : if(sscanf(line,"%s%lf%lf%lf",name,&a.x,&a.y,&a.z)<4) return(0);
                    Get_Type(name,type);
                    break;
  }
  name[ATM_LEN]=0;
  type[ATM_LEN]=0;

  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,name);
  M->atoms->fv=VequalV(a);

  return(1);
}

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

int Get_Symmetry_Operator(Mol *M, char *line, unsigned mode, unsigned type)
{
  int i;
  Matrix T={{1.0,0.0,0.0},{0.0,1.0,0.0},{0.0,0.0,1.0}};
  Vector t={0.0,0.0,0.0};
  Symm *new=NULL;

  if(M==NULL) return(0);
  if(line==NULL) return(0);
  if(strlen(line)>=SYM_LEN) line[SYM_LEN-1]=0;
  if(strchr(line,'\n')!=NULL) *(strchr(line,'\n'))=0;
  i=strlen(line);
  while(i--,line[i]==' ') line[i]=0;
  Last_Character(line,' ');
  if(mode==F_GSTCOR)
  {
    if(!Get_Symmetry_Matrix(line,&T,&t)) return(0);
  }
  else if(!Get_Symmetry(line,&T,&t)) return(0);

  if((new=New_Symm(M->symmetry))==NULL) return(0);
  M->symmetry=new;
  Initialize_Symmetry(M->symmetry,line,T,t,type);

  if(!Identity_Operator(M->symmetry))
  {
    if(type==S_NORM) M->geneq++;
    else if(type==S_CENT) M->geneq*=2;
  }
  return(1);
}

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

int Shelx_Card_Type(char *line)
{
  if(!strncmp(line,"ZERR",4))  return(BAD);
  if(!strncmp(line,"SFAC",4))  return(BAD);
  if(!strncmp(line,"UNIT",4))  return(BAD);
  if(!strncmp(line,"L.S.",4))  return(BAD);
  if(!strncmp(line,"ACTA",4))  return(BAD);
  if(!strncmp(line,"FREE",4))  return(BAD);
  if(!strncmp(line,"FMAP",4))  return(BAD);
  if(!strncmp(line,"GRID",4))  return(BAD);
  if(!strncmp(line,"FVAR",4))  return(BAD);
  if(!strncmp(line,"HKLF",4))  return(BAD);
  if(!strncmp(line,"AFIX",4))  return(BAD);
  if(!strncmp(line,"DFIX",4))  return(BAD);
  if(!strncmp(line,"TITL",4))  return(TIT);
  if(!strncmp(line,"CELL",4))  return(CEL);
  if(!strncmp(line,"LATT",4))  return(LAT);
  if(!strncmp(line,"SYMM",4))  return(SYM);
  if(!strncmp(line,"END",3))   return(END);
  return(ATM);
}

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

int End_Fractional_Coordinates(Mol *M)
{
  Atms *A=NULL;
  if(M==NULL) return(0);
  if(!(M->mode&FRAC_BIT)) return(1);
  Error_Message(E_ENDFRC,"End Fractional Coordinates");
  M->mode&=(FULL_BIT^FRAC_BIT);
  M->symmetry=Close_All_Symmetries(M->symmetry);
  for(A=First_Atom(M->atoms);A!=NULL;A=A->next) A->fv=VequalV(A->v);
  M->a=1.0;
  M->b=1.0;
  M->c=1.0;
  M->al=90.0;
  M->be=90.0;
  M->ga=90.0;
  M->cellvol=0.0;
  M->freevol=0.0;
  M->geneq=1;
  M->centring='p';
  return(1);
}

/***************************************************************************/
/*******************  Crystal Calculations  ********************************/
/***************************************************************************/

int Calculate_Free_Unit_Cell_Volume(Mol *M, char *line, Set *set)
{
  Vector min={0.0,0.0,0.0}, max={0.0,0.0,0.0}, d={0.0,0.0,0.0};
  double V;
  long total;
  double tot_val=0.0;
  double err_val=0.0;
  int gnum=0;

  Box_Atoms(M,"1 1 1",set);

/* define encompassing box */

  Find_Encompassing_Box(M,FULL_BIT,&min,&max,&d);

/* choose two ways of calulating volume */

  if(set->mode&MTC_BIT)   /* Monte Carlo approach */
  {
    tot_val=Monte_Carlo_Volume(M,set,gnum,&min,&d,VOL_UC);
    err_val=set->veps;
  }
  else                    /* systemmatic approach */
  {
    tot_val=Fixed_Grid_Volume(M,set,gnum,&min,&max,&d,&total,VOL_UC);
    err_val=V/(double)total;
  }
  M->freevol=M->cellvol-tot_val;
  sprintf(line,"Free Unit Cell Volume =%10.2f cubic Angstroms [Cell=%10.2f Atoms=%10.2f]"
              ,M->freevol,M->cellvol,tot_val);
  Out_Message(line,O_NEWLN);
  return(1);
}

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

int Calculate_Group_Volume(Mol *M, char *argline, Set *set, unsigned mode)
{
  int i,n;
  double value=0.0;
  char *args[MAXARG];
  int num=0;
  Grps *group;
  char line[LINELEN];
  char outline[LINELEN];

  if(M->groups==NULL)
  {
    Error_Message(E_NOGRPS,"Calculate Group Volume");
    return(0);
  }

/* load specified group names into memory                                 */

  if(strlen(argline)<1) return(0);
  strcpy(line,argline);
  if((num=Get_Group_Names(line,outline,args))<1) return(0);

  if((mode==CAVV)&&(set->modev&ABOX_BIT)) Box_Atoms(M,"-1",set);

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

  for(n=0;ngroups);group!=NULL;group=group->next,i++)
    {
      if(stricmp(group->name,args[n])==0)
      {
        if(mode==CAVV)
        {
          sprintf(line,"Calculating Volume of cavity of group %d (%s) at R=%4.1f"
                      ,i,group->name,group->volrad);
          Out_Message(line,O_NEWLN);

          group->mode|=CAVV_BIT;
          value=Molecular_Volume(M,NULL,set,group);
          group->mode&=(FULL_BIT^CAVV_BIT);

          sprintf(line,"%s Group %s R=%4.1f cavity volume = %f cubic angstroms"
                      ,M->name,group->name,group->volrad,value);
          Out_Message(line,O_NEWLN);
        }
        else
        {
          sprintf(line,"Calculating Volume of group %d (%s)",i,group->name);
          Out_Message(line,O_NEWLN);
          value=Molecular_Volume(M,NULL,set,group);
          sprintf(line,"%s Group %s volume = %f cubic angstroms"
                      ,M->name,group->name,value);
          Out_Message(line,O_NEWLN);
        }
        Out_Message("",O_NEWLN);
      }
    }
  }
  return(1);
}

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

int Exclude_Atoms_Outside_Groups(Mol *M, Set *set, char *argline)
{
  Atms *at=NULL;
  int i,n,g;
  char *args[MAXARG];
  int gnum[MAXARG];
  Grps *G[MAXARG];
  double volrad[MAXARG];
  double atsep;
  int num=0;
  Grps *group=NULL;
  char line[LINELEN];
  char outline[LINELEN];

  if(M->groups==NULL)
  {
    Error_Message(E_NOGRPS,"Calculate Group Volume");
    return(0);
  }

/* load specified group names into memory                                 */

  if(strlen(argline)<1) return(0);
  strcpy(line,argline);
  if((num=Get_Group_Names(line,outline,args))<1) return(0);

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

  sprintf(line,"Closing atoms outside groups");

  for(g=0,n=0;ngroups);group!=NULL;group=group->next,i++)
    {
      if(stricmp(group->name,args[n])==0)
      {
        if(strlen(line)volrad;
        g++;
        if(g>=MAXARG-1) break;
      }
    }
  }
  
  if(g) Out_Message(line,O_NEWLN);
  else
  {
    Out_Message("No matching groups found",O_NEWLN);
    return(0);
  }

/* delete all atoms not within specified groups range                      */

  for(at=First_Atom(M->atoms);at!=NULL;at=at->next)
  {
    at->bond=Close_All_Bonds(at->bond);
    at->bonds=0;
  }

  at=First_Atom(M->atoms);
  while((g!=0)&&(at!=NULL))
  {
    if(at==NULL) break;
    else M->atoms=at;
    if(strcmp(at->type,"A")==0) {at=at->next; continue; }
    for(i=0;imain==NULL)) continue;
      atsep=Vmag(Vdif(at->v,G[i]->main->v))-at->radius-volrad[i];
      if(atsep<=0.0) break;
    }
    if(i==g) at=Close_Current_Atom(at);
    else at=at->next;
  }

  for(i=0,at=First_Atom(M->atoms);at!=NULL;at=at->next,i++);
  M->num_atoms=i;

/* refind all bonded groups */

  for(i=0;i