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

/**************************************************************************/
/******************       Group Functions        **************************/
/**************************************************************************/
/******************        This module is        **************************/
/******************      system independant      **************************/
/**************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "sterdefn.h" /* declaration of structures, functions and globals */
#include "stercomm.h" /* definitions of all menu command options          */
#include "sterfile.h"      /* functions for all file manipulations        */
#include "stermem.h"  /* dynamic memory management module                 */
#include "stertext.h" /* all text functions for text mode                 */
#include "steraid.h"  /* additional functions needed                      */
#include "stergrp.h"  /* group manipulation functions                     */

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

void Initialize_Group(Grps *group)
{
  group->name[0]=0;         /* group lable for identification             */
  group->atex[0]=0;         /* list of extra atom types allowed           */
  group->main=NULL;         /* atom no of atom bonded to origin           */
  group->num_atoms=0;       /* no of atoms in group                       */
  group->SVangle=0.0;       /* peak semi vertex angle for group           */
  group->volume=0.0;        /* group volume                               */
  group->cavity=0.0;        /* cavity volume                              */
  group->radvol=0.0;        /* radial profile volume                      */
  group->radcav=0.0;        /* radvol cavity volume                       */
  group->volrad=0.0;        /* radius for cavity volume calculation       */
  group->mode=FIXA_BIT;     /* mode (fixed numatoms or not)               */
/* next and prev pointers have already been assigned, do not fiddle !!!   */
}

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

int Get_Group_Names(char *inline, char *outline, char **oargs)
{
  int i,n,a;
  char *iargs[MAXARG];
  char *args[MAXARG];
  int inum=0,onum=0,num=0,count=0;
  FILE *GR=NULL;
  char line[LINELEN];
  char name[LINELEN];

  sprintf(line,"steric.grp");
  if((GR=fopen(line,"rt"))==NULL)
  {
    sprintf(line,"%s/steric.grp",STERICHOME);
    if((GR=fopen(line,"rt"))==NULL)
    {
      sprintf(line,"Get Group Names [steric.grp] | [%s/steric.grp]",STERICHOME);
      Error_Message(E_NOGRPF,line);
      return(0);
    }
  }
  if((strlen(inline)<1)||((inum=Get_Arguements(inline,iargs))<1)) return(0);
  outline[0]=0;

/* load specified group names into memory                                 */

  for(a=0,i=0,n=0;n<inum;i=0,n++)
  {
    rewind(GR);
    strcpy(name,iargs[n]);
    while((get_next_line(GR,"GROUP",line)>0L)&&(!feof(GR)))
    {
      if(strlen(outline)>LINELEN-10) break;
      if(strchr(line,'\n')!=NULL) strchr(line,'\n')[0]=0;
      if((num=Get_Arguements(line,args))>0)
      for(count=1;count<num;count++)
      {
        if(stricmp(name,args[count])==0)
        {
          strcpy(name,args[num-1]);
          i++; a++;
          strcat(outline,name);
          strcat(outline," ");
          break;
        }
      }
      if(i) break;
    }
    if(i==0)
    {
      i++; a++;
      strcat(outline,name);
      strcat(outline," ");
      break;
    }
    if(a>=MAXARG) break;
  }
  fclose(GR);
  onum=Get_Arguements(outline,oargs);
  return(onum);
}

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

int Get_Groups_by_Name(Mol *M, char **args, int num, Grps **group)
{
  int i,n;
  int gi=0;

  if(M==NULL) return(0);
  if(M->groups==NULL) return(0);
  for(i=0,n=0;n<num;n++)
  {
    if(i>=MAXARG) break;
    group[i]=NULL;
    if(((group[i]=Goto_Group_Name(M->groups,args[n]))!=NULL)
     ||((gi=Int_String(args[n]))!=0))
    {
      if((group[i]==NULL)&&(gi)) group[i]=Goto_Group(M->groups,gi);
    }
    if(group[i]!=NULL) i++;
    else
    {
      Error_Message(E_NOFGRP,"Get_Groups_by_Name");
      continue;
    }
  }
  return(i);
}

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

int Same_Bonds(Atms *A, Atms *B, int dif)
{
  Bond *bdA=NULL, *bdB=NULL;
  int a,b;
  if(dif>1000) return(1);
  for(b=0,bdB=First_Bond(B->bond);bdB!=NULL;bdB=bdB->next,b++);
  for(a=0,bdA=First_Bond(A->bond);bdA!=NULL;bdA=bdA->next,a++);
  if((a-b!=dif)&&(a-b!=0)) return(0);
  for(a=0,b=0,bdB=First_Bond(B->bond);bdB!=NULL;bdB=bdB->next,b++)
  {
    for(bdA=First_Bond(A->bond);bdA!=NULL;bdA=bdA->next)
    {
      if(stricmp(bdA->atom->type,bdB->atom->type)==0)
      {
        a++;
        break;
      }
    }
  }
  if(a==b) return(1);
  return(0);
}

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

int Same_Atom(Atms *A, Atms *B, int dif)
{
  if(stricmp(A->type,B->type)!=0) return(0);
  if(!Same_Bonds(A,B,dif)) return(0);
  return(1);
}

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

int Same_Types(Grps *mgroup, Grps *group)
{
  Atms *at=NULL;
  char *args[MAXARG];
  char line[LINELEN];
  int num=0,i,n;
  
  if((mgroup==NULL)||(group==NULL)) return(0);
  if(group->mode&FIXA_BIT) return(0);
  if(mgroup->num_atoms<group->num_atoms) return(0);
  n=Get_Group_Number(mgroup);
  sprintf(line,"%s ",group->atex);
  for(at=First_Atom(group->main);at!=NULL;at=at->next)
  {
    strcat(line,at->type);
    if(strlen(line)>LINELEN-5) break;
    strcat(line," ");
  }
  if((num=Get_Arguements(line,args))<1) return(0);

  for(at=First_Atom(mgroup->main);at!=NULL;at=at->next)
  {
    if(at->group!=n) continue;
    for(i=0;i<num;i++)
    {
      if(stricmp(at->type,args[i])==0) break;
    }
    if(i==num) return(0);
  }
  return(1);
}

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

int Group_Matches(Grps *mgroup, Grps *group)
{
  Atms *atM=NULL, *atG=NULL;
  int a,b,d, mgn;
  
  if((group==NULL)||(mgroup==NULL)) return(0);
  if(stricmp(mgroup->main->type,group->main->type)!=0) return(0);
  if((Same_Types(mgroup,group))||(group->num_atoms==mgroup->num_atoms))
  {
    for(a=0,b=0,d=1,atG=First_Atom(group->main);atG!=NULL;atG=atG->next,a++,d=0)
    {
      if(!(group->mode&FIXA_BIT)) d=1001;
      for(atM=First_Atom(mgroup->main);atM!=NULL;atM=atM->next)
      {
        mgn=Get_Group_Number(mgroup);
        if((atM->group==mgn)&&(Same_Atom(atM,atG,d)))
        {
          b++;
          break;
        }
      }
    }
    if(a==b) return(1);
  }
  return(0);
}

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

Grps *Load_Group(FILE *GR, Grps *group, char *line, char *name)
{
  char *args[MAXARG];
  int num=0,i;
  Atms *at=NULL;
  int a, b;
  Atms *atomA=NULL, *atomB=NULL;

  if((GR==NULL)||(strncmp(line,"GROUP",5)!=0)) return(NULL);

  group=New_Group(group);
  Initialize_Group(group);
  strcpy(group->name,name);
  while(fgets(line,148,GR),!feof(GR))
  {
    if(line[0]=='#') break;
    if(strchr(line,'\n')!=NULL) strchr(line,'\n')[0]=0;
    if((num=Get_Arguements(line,args))>0)
    if(stricmp(args[0],"GROUP")==0) break;
    if(stricmp(args[0],"ATOM")==0)
    {
      for(i=1;i<num;i++)
      {
        at=New_Atom(at);
        Initialize_Atom(at);
        strncpy(at->name,args[i],ATM_LEN);
        strncpy(at->type,args[i],ATM_LEN);
        group->num_atoms++;
      }
      group->main=First_Atom(at);
    }
    else if(stricmp(args[0],"ATEX")==0)
    {
      for(i=1;i<num;i++)
      {
        if(i>1) strcat(group->atex," ");
        strcat(group->atex,args[i]);
      }
      group->mode&=FULL_BIT^FIXA_BIT;
    }
    else if((stricmp(args[0],"VOLR")==0)&&(num>1))
    {
      sscanf(args[1],"%lf",&group->volrad);
      if(group->volrad<0.0) group->volrad=0.0;
      if(group->volrad>MAXVRAD) group->volrad=MAXVRAD;
    }
    else if((stricmp(args[0],"BOND")==0)&&(num>2))
    {
      a=0;b=0;
      sscanf(args[1],"%d",&a);
      atomA=Goto_Atom(group->main,a);
      for(i=2;i<num;i++)
      {
        sscanf(args[i],"%d",&b);
        atomB=Goto_Atom(group->main,b);
        Add_New_Bond(atomA,atomB,SINGLE_B);
      }
    }
  }
  if(group->volrad<0.0) group->volrad=0.0;
  return(group);
}

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

int Kill_All_Groups(Mol *M, char *name)
{
  char *args[MAXARG];
  int num=0,count=0,i=0;
  Atms *at=NULL;
  Grps *mgroup=NULL;
  FILE *GR=NULL;
  char line[LINELEN];
  int a;

  if(M->groups==NULL)
  {
    Error_Message(E_NOGRPS,"Kill All Groups");
    return(0);
  }
  if((GR=fopen("steric.grp","rt"))==NULL)
  {
    Error_Message(E_NOGRPF,"Kill All Groups");
    return(0);
  }

/* load specified group structure into memory                             */

  while(get_next_line(GR,"GROUP",line)>0L)
  {
    if(feof(GR)) return(0);
    if(strchr(line,'\n')!=NULL) strchr(line,'\n')[0]=0;
    if((num=Get_Arguements(line,args))>0)
    for(count=1;count<num;count++)
    {
      if(stricmp(name,args[count])==0)
      {
        strcpy(name,args[num-1]);
        i++;
        break;
      }
    }
  }
  fclose(GR);
  if(i==0)
  {
    Error_Message(E_UNKGRP,"Kill All Groups");
  }
  if(i>1)
  {
    Error_Message(E_TOMGRP,"Kill All Groups");
  }

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

  a=0;
  for(mgroup=First_Group(M->groups);mgroup!=NULL;mgroup=mgroup->next)
  {
    if(stricmp(mgroup->name,name)==0)
    {
      a++;
      for(at=First_Atom(M->atoms);at!=NULL;at=at->next)
      {
        if(at->group==Get_Group_Number(mgroup))
          at->stat=at->stat&(!MAIN_BIT);
      }
    }
  }
  if(a)
  {
    sprintf(line,"Killed %d %s groups",a,name);
    Out_Message(line,O_NEWLN);
  }

  return(1);
}

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

int Name_All_Groups(Mol *M)
{
  char *args[MAXARG];
  char name[LINELEN];
  int num=0;
  Grps *new=NULL, *group=NULL,*mgroup=NULL, *lgroup=NULL;
  FILE *GR=NULL;
  char line[LINELEN];

  if(M->groups==NULL)
  {
    Error_Message(E_NOGRPS,"Name All Groups");
    return(0);
  }
  if((GR=fopen("steric.grp","rt"))==NULL)
  {
    Error_Message(E_NOGRPF,"Name All Groups");
    return(0);
  }

/* load all group structures into memory                                  */

  while(get_next_line(GR,"GROUP",line)>0L)
  {
    if(feof(GR)) return(0);
    if(strchr(line,'\n')!=NULL) strchr(line,'\n')[0]=0;
    if((num=Get_Arguements(line,args))>0)
    strcpy(name,args[num-1]);
    if((new=Load_Group(GR,lgroup,line,name))==NULL)
    {
      Error_Message(E_LDGRP,"Name All Groups");
      break;
    }
    lgroup=new;
  }
  fclose(GR);

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

  for(group=First_Group(lgroup);group!=NULL;group=group->next)
  {
    for(mgroup=First_Group(M->groups);mgroup!=NULL;mgroup=mgroup->next)
    {
      if(Group_Matches(mgroup,group))
      {
        sprintf(line,"Renaming group %s to %s",mgroup->name,group->name);
        Out_Message(line,O_NEWLN);
        strcpy(mgroup->name,group->name);
        if(!AlmostZero(group->volrad)) mgroup->volrad=group->volrad;
      }
    }
  }

/* close down temporary group memory                                      */

  for(group=First_Group(lgroup);group!=NULL;)
  {
    Close_All_Atoms(group->main,0);
    group=Close_Group(group);
  }

  return(1);
}

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

double Find_Max_Group_SVangles(Mol *M)
{
  double cone=0.0;
  Atms *at=NULL;
  Grps *group=NULL;
  int n=0,p=0,imax=0;

  if(M->groups==NULL)
  {
    for(n=0,at=First_Atom(M->atoms);at!=NULL;at=at->next)
    {
      n++; if(n>M->num_atoms) Error_Message(E_ATMNUM,"Find Max Group SVangles");
      if((at->stat&MAIN_BIT)&&(!AlmostZero(at->SVangle))&&(cone<at->theta+at->SVangle))
      {
        imax=n;
        cone=at->theta+at->SVangle;
      }
    }
    p=1;
  }
  else
  {
    group=First_Group(M->groups);
    if(group->name[0]=='0')
    {
      group->SVangle=group->main->SVangle;
      group=group->next;
    }
    for(p=0;group!=NULL;group=group->next,p++)
    {
      group->SVangle=0.0;
      for(n=0,at=First_Atom(M->atoms);at!=NULL;at=at->next)
      {
        n++; if(n>M->num_atoms) Error_Message(E_ATMNUM,"Find Max Group SVangles");
        if((at->stat&MAIN_BIT)&&((at->group==p+1)||(at->group==0))
         &&(!AlmostZero(at->SVangle))&&(group->SVangle<at->theta+at->SVangle))
        {
          group->SVangle=at->theta+at->SVangle;
          imax=n;
        }
      }
      cone+=group->SVangle;
    }
    if(p<3)
    {
      if(group=First_Group(M->groups),group->name[0]=='0')
      {
        p++;
        cone+=group->SVangle;
      }
    }
  }
  M->atoms=Goto_Atom(M->atoms,imax);
  return(cone*2/p);
}

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

int Find_Bonded_Atoms(Grps *group, Atms *o, int i)
{
  Bond *bond=NULL;
  
  if((o!=NULL)&&(!(o->stat&GRP_BIT))&&(o->stat&MAIN_BIT))
  {
    o->stat|=GRP_BIT;
    o->group=i;
    group->num_atoms++;
    if(group->main==NULL) group->main=o;
    for(bond=First_Bond(o->bond);bond!=NULL;bond=bond->next)
    {
      Find_Bonded_Atoms(group,bond->atom,i);
    }
    return(1);
  }
  return(0);
}

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

int Assign_Groups(Mol *M)
{
  int i;
  Bond *bond=NULL;
  Grps *new=NULL;
  Atms *atom=NULL;
  Atms *o=M->origin;
  
  M->groups=Close_All_Groups(M->groups);
  if(o==NULL) 
  {
    M->groups=Close_All_Groups(M->groups);
    return(0);
  }
  for(atom=First_Atom(o);atom!=NULL;atom=atom->next)
  {
    atom->stat&=255^GRP_BIT;
    atom->group=0;
  }
  o->stat|=GRP_BIT;
  for(i=0,bond=First_Bond(o->bond);bond!=NULL;bond=bond->next)
  {
    if((bond->atom!=NULL)&&(bond->atom->stat&MAIN_BIT))
    {
      i++;
      atom=bond->atom;
    }
  }
  if(i<1)
  {
    M->groups=Close_All_Groups(M->groups);
    return(0);
  }
  else if(i==1)
  {
    if(atom!=NULL)
    {
      o=atom;
      o->stat|=GRP_BIT;
      M->main_atom=o;
      if((new=New_Group(M->groups))==NULL) return(0);
      M->groups=new;
      Initialize_Group(M->groups);
      M->groups->name[0]='0';
      M->groups->name[1]=0;
      M->groups->num_atoms=1;
      M->groups->SVangle=o->SVangle;
      M->groups->main=o;
    }
  }
  for(i=1,bond=First_Bond(o->bond);bond!=NULL;bond=bond->next)
  {
    if((new=New_Group(M->groups))==NULL) return(0);
    M->groups=new;
    Initialize_Group(M->groups);
    sprintf(M->groups->name,"%d",i);
    if(!Find_Bonded_Atoms(M->groups,bond->atom,i))
      M->groups=Close_Group(M->groups);
    else i++;
  }
  Name_All_Groups(M);
  return(1);
}

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

int Find_All_Bonded_Groups(Mol *M, unsigned mode)
{
  int i;
  Grps *new=NULL;
  Atms *atom=NULL;
  char line[LINELEN];
  
  M->groups=Close_All_Groups(M->groups);
  for(atom=First_Atom(M->atoms);atom!=NULL;atom=atom->next)
  {
    if(atom->group==FULL_BIT)
    {
      atom->stat|=GRP_BIT;
      continue;
    }
    if((mode==NO_SINGLE)&&(atom->bonds==0))
    {
      atom->stat|=GRP_BIT;
      continue;
    }
    atom->stat&=FULL_BIT^GRP_BIT;
    atom->group=0;
  }
  for(i=1,atom=First_Atom(M->atoms);atom!=NULL;atom=atom->next)
  {
    if(atom->stat&GRP_BIT) continue;
    if((new=New_Group(M->groups))==NULL) return(0);
    M->groups=new;
    Initialize_Group(M->groups);
    sprintf(M->groups->name,"%d",i);
    if(!Find_Bonded_Atoms(M->groups,atom,i))
      M->groups=Close_Group(M->groups);
    else
    {
      sprintf(line,"Found group %s of %d atoms",M->groups->name,M->groups->num_atoms);
      Out_Message(line,O_NEWLN);
      i++;
    }
  }
  Name_All_Groups(M);
  return(1);
}

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