#include "utility.h"
#include "newgeo.h"

int intgen()
/*============================================================================*/
/* PURPOSE: GENERATE THE INTERNAL COORDINATES FOR THE CURRENT GEOMETRY.
*/
   {
   if (!continuously_bonded()) return(FALSE);
   else if (!generate_reference_atoms()) return(FALSE);
   else if (!generate_internal_coords()) return(FALSE);
   else return(TRUE);
   }

int continuously_bonded()
/*============================================================================*/
/*  PURPOSE: VERIFY THAT ALL ATOMS ARE CONTINUOUSLY BONDED.  THIS IS NECESSARY
/*	ONLY BECAUSE THE INTERNAL COORDINATE GENERATION ROUTINE FROM JEFF GAW
/*	WILL FATALLY FAIL ON NON CONNECTED MOLECULES.
*/
   {
   int atom_reachable[MAX_ATOMS];	/* T/F: ATOM REACHABLE FROM ROOT. */
   int natom;				/* WHICH ATOM. */

   void mark_connected_atoms();		/* MARK ALL CONNECTED ATOMS. */

   for (natom=0; natom<atom_cnt; natom++)
      {
      atom_reachable[natom]=FALSE;	/* NOTHING BONDED TO START. */
      }

   /* MARK ALL ATOMS THAT ARE CONNECTED TO ATOM 0. */
   mark_connected_atoms(0,atom_reachable);

   /* VERIFY THAT ALL ATOMS WERE CONNECTED TO ATOM 0. */
   for (natom=0; natom<atom_cnt; natom++)
      {
      if (!atom_reachable[natom]) return(errmsg("This molecule is not completely connected",FALSE));
      }

   return(TRUE);	/* MOLECULE MUST BE COMPLETELY CONNECTED. */
   }

void mark_connected_atoms(this_atom,atom_bonded)
/*============================================================================*/
/* PURPOSE: TO MARK ALL ATOMS CONNECTED TO THIS_ATOM.
*/
   int this_atom;
   int atom_bonded[MAX_ATOMS];
   {
   int nbond;	/* WHICH BOND. */

   if (!atom_bonded[this_atom])
      {
      atom_bonded[this_atom]=TRUE;	/* MARK THIS ATOM. */
      for (nbond=0; nbond<intco[this_atom].bond_cnt; nbond++)	/* MARK ALL CONNECTED ATOMS. */
         {
         mark_connected_atoms(intco[this_atom].bonds[nbond],atom_bonded);
         }
      }

   return;
   }

int generate_reference_atoms()
/*============================================================================*/
   {
   int i;
   int j;
   int max_con;
   int parent;
   int gparent;
   int ggparent;
   int assign_ggparent;
   int new_order[MAX_ATOMS];
   int assigned;
   int look_at;
   char reorder_file_name[130];		/* NAME OF ATOM REORDER FILE. */
   FILE *reorder_file;
   INTERNAL_COORDINATES swap[MAX_ATOMS];

   /* REORDER THE ATOMS BASED ON THE NUMBER OF CONNECTIONS */
   max_con = 0;
   for (i=0; i<atom_cnt; i++)
      {
      if (intco[i].bond_cnt > max_con)
         {
         j = i;
         max_con = intco[i].bond_cnt;
         }
      }
   new_order[0] = j;

   assigned = 1;
   look_at = 0;
   while (look_at < atom_cnt)
      {
      while((j=next_bonded(new_order[look_at],new_order,assigned)) >= 0)
         {
         new_order[assigned++] = j;
         }
      look_at++;
      }

   /* TELL USER THE NEW ORDER OF THINGS SO THAT THEY CAN RE-ARRANGE THE WORLD TO CORRESPOND. */
   strcpy(reorder_file_name,input_file_name);
   new_extension(reorder_file_name,".reorder");
   reorder_file = fopen(reorder_file_name,"w");
   for (i=0; i<atom_cnt; i++)
      {
      fprintf(reorder_file,"Old atom number %2d (%.4s) is now number %2d\n",
         new_order[i]+1,intco[new_order[i]].atom_name,i+1);
      }
   fclose(reorder_file);

   /* NOW, PUT THINGS IN THE CORRECT ORDER */
   for (i=0; i<atom_cnt; i++)
      {
      memcpy(&swap[i],&intco[new_order[i]],sizeof(INTERNAL_COORDINATES));
      }
   memcpy(intco,swap,atom_cnt*sizeof(INTERNAL_COORDINATES));

   /* FIX THE COORDINATES SO THAT THE FIRST ATOM IS AT THE ORIGIN */
/*   for (i=1; i<atom_cnt; i++)
      {
      for (j=0; j<3; j++)
         {
         intco[i].vector[j] -= intco[0].vector[j];
         }
      }
   intco[0].vector[0] = 0.0;
   intco[0].vector[1] = 0.0;
   intco[0].vector[2] = 0.0;
*/

   /* FIX THE CONNECTION NUMBERS TO MATCH THE NEW ATOM ORDER */
   for (i=0; i<atom_cnt; i++)
      {
      for (j=0; j<intco[i].bond_cnt; j++)
         {
         for (look_at=0; look_at<atom_cnt; look_at++)
            {
            if (new_order[look_at] == intco[i].bonds[j]) break;
            }
         intco[i].bonds[j] = look_at;
         }
      }

   /* FIX THE COORDINATES OF THE Z_SLOT TO BE THE UNIT VECTOR TO NEGATIVE */
   /* Z AXIS (USED FOR THE GPARENT AND GGPARENT COORDINATES */
   strncpy(intco[Z_SLOT].atom_name,"Z    ",5);
/*   intco[Z_SLOT].vector[0] = -intco[1].vector[0];
   intco[Z_SLOT].vector[1] = -intco[1].vector[1];
   intco[Z_SLOT].vector[2] = -intco[1].vector[2];*/
   intco[Z_SLOT].vector[0] = 2.0 * intco[0].vector[0] - intco[1].vector[0];
   intco[Z_SLOT].vector[1] = 2.0 * intco[0].vector[1] - intco[1].vector[1];
   intco[Z_SLOT].vector[2] = 2.0 * intco[0].vector[2] - intco[1].vector[2];

   assign_ggparent = X_SLOT;
/*   if (intco[Z_SLOT].vector[1]==0 && intco[Z_SLOT].vector[2]==0) assign_ggparent=Y_SLOT;*/
   strncpy(intco[X_SLOT].atom_name,"X    ",5);
/*   intco[X_SLOT].vector[0] = 1.0 - intco[0].vector[0];
   intco[X_SLOT].vector[1] = 0.0 - intco[0].vector[1];
   intco[X_SLOT].vector[2] = 0.0 - intco[0].vector[2];*/
   intco[X_SLOT].vector[0] = intco[0].vector[0] + 1.0;
   intco[X_SLOT].vector[1] = intco[0].vector[1] + 0.0;
   intco[X_SLOT].vector[2] = intco[0].vector[2] + 0.0;
   if (intco[Z_SLOT].vector[1]==intco[X_SLOT].vector[1] && intco[Z_SLOT].vector[2]==intco[X_SLOT].vector[2]) assign_ggparent=Y_SLOT;

   strncpy(intco[Y_SLOT].atom_name,"Y    ",5);
/*   intco[Y_SLOT].vector[0] = 0.0;
   intco[Y_SLOT].vector[1] = 1.0;
   intco[Y_SLOT].vector[2] = 0.0;*/
   intco[Y_SLOT].vector[0] = intco[0].vector[0] + 0.0;
   intco[Y_SLOT].vector[1] = intco[0].vector[1] + 1.0;
   intco[Y_SLOT].vector[2] = intco[0].vector[2] + 0.0;

   /* DETERMINE THE REFERENCE ATOMS. */
   intco[0].ref[0] = intco[0].ref[1] = intco[0].ref[2] = -1;

   intco[1].ref[0] = 0;
   intco[1].ref[1] = Z_SLOT;
   intco[1].ref[2] = -1;

   for (i=2; i<atom_cnt; i++)
      {
      /* SELECT A PARENT FOR THE CURRENT ATOM. */
      parent = X_SLOT;
      for (j=0; j<intco[i].bond_cnt; j++)
         {
         if (intco[i].bonds[j]<parent) parent=intco[i].bonds[j];
         }

      /* SELECT A GRANDPARENT FOR THE CURRENT ATOM. */
      gparent = X_SLOT;
      for (j=0; j<intco[parent].bond_cnt; j++)
         {
         if (intco[parent].bonds[j]<gparent && intco[parent].bonds[j]<parent) gparent=intco[parent].bonds[j];
         }

      /* SELECT A GREAT-GRANDPARENT FOR THE CURRENT ATOM. */
      if (gparent < Z_SLOT)
         {
         ggparent = assign_ggparent;
         for (j=0; j<intco[gparent].bond_cnt; j++)
            {
            if (intco[gparent].bonds[j]<ggparent && intco[gparent].bonds[j]<gparent)
               ggparent=intco[gparent].bonds[j];
            }

         /* IF THAT DID NOT FIND ANYTHING, THEN LOOK FOR SOMETHING WHICH IS */
         /* CONNECTED TO GPARENT WHICH IS NOT PARENT AND LOWER NUMBER */
         if (ggparent == assign_ggparent)
            {
            for (j=0; j<intco[gparent].bond_cnt; j++)
               {
               if (intco[gparent].bonds[j]!=parent && intco[gparent].bonds[j]<i)
                  {
                  ggparent = intco[gparent].bonds[j];
                  break;
                  }
               }
            }
         }
      else
         {
         gparent = 1;
         if (i == 2) ggparent = assign_ggparent;
         else
            {
            ggparent = assign_ggparent;
            for (j=0; j<intco[parent].bond_cnt; j++)
               {
               if (intco[parent].bonds[j]!=gparent && intco[parent].bonds[j]<i)
                  {
                  ggparent = intco[parent].bonds[j];
                  break;
                  }
               }
            }
         }
/*DEBUG
      printf("intco[%d].ref[0]=%d, .ref[1]=%d, .ref[2]=%d\n",i,parent,gparent,ggparent);
/**/
      intco[i].ref[0] = parent;
      intco[i].ref[1] = gparent;
      intco[i].ref[2] = ggparent;
      }

   return(TRUE);
   }

int next_bonded(from,new_order,assigned)
/*============================================================================*/
/* PURPOSE: SELECT THE MAXIMALLY BONDED ATOM CONNECTED TO "FROM" NOT ALREADY
/*	IN THE "ASSIGNED" ELEMENTS OF "NEW_ORDER".
*/
    int from;
    int new_order[MAX_ATOMS];
    int assigned;
    {
    int max_con;
    int i,j,k,found;
    int check_atom;

    max_con = 0;
    j = -1;

    for (i=0; i<intco[from].bond_cnt; i++)
	{
        check_atom=intco[from].bonds[i];

        found=FALSE;
	for (k=0; k<assigned && !found; k++)
	    {
	    if (new_order[k] == check_atom) found = TRUE;
	    }

	/* IF IT IS NOT IN THE ASSIGNED LIST, SEE IF IT HAS MORE CONNECTIONS */
	if (!found && (k=intco[check_atom].bond_cnt) > max_con)
	    {
	    max_con = k;
	    j = intco[from].bonds[i];
	    }
	}

    return(j);
    }

int generate_internal_coords()
/*============================================================================*/
   {
   int i;
   int parent;
   int gparent;
   int ggparent;
   double bond1[3];
   double bond2[3];
   double bond3[3];
   double bond4[3];
   double cross1[3];
   double cross2[3];
   double cross3[3];
   double length1;
   double length2;

   /* CALCULATE ALL THE BOND LENGTHS (WHICH START AT ATOM ONE). */
   for (i=1; i<atom_cnt; i++)
      {
      parent = intco[i].ref[0];
      vec_diff(intco[i].vector,intco[parent].vector,bond1);
      intco[i].length=vec_length(bond1);
      } 

   /* NOW THE ANGLES (WHICH START AT ATOM TWO). */
   for (i=2; i<atom_cnt; i++)
      {
      parent = intco[i].ref[0];
      gparent = intco[i].ref[1];
      vec_diff(intco[parent].vector,intco[i].vector,bond1);
      vec_diff(intco[parent].vector,intco[gparent].vector,bond2);
      intco[i].angle=vec_angle(bond1,bond2);
      }

   /* NEXT COME THE DIHEDRALS (WHICH START AT ATOM TWO). */
   for (i=2; i<atom_cnt; i++)
      {
      parent = intco[i].ref[0];
      gparent = intco[i].ref[1];
      ggparent = intco[i].ref[2];
      vec_diff(intco[parent].vector,intco[i].vector,bond1);
      vec_diff(intco[parent].vector,intco[gparent].vector,bond2);
      vec_diff(intco[gparent].vector,intco[parent].vector,bond3);
      vec_diff(intco[gparent].vector,intco[ggparent].vector,bond4);
      vec_crossprod(bond1,bond2,cross1);
      vec_crossprod(bond3,bond4,cross2);

      /* CHECK THE DEGENERATE CASE OF NO TORSION ANGLE - LINEAR CONNECTION. */
      length1=vec_length(cross1);
      length2=vec_length(cross2);
      if (length1==0 || length2==0)
         {
/*DEBUG
         printf("forced LC;length1=%f,length2=%f,atom=%d,parent=%d,gparent=%d,ggparent=%d\n",length1,length2,i,parent,gparent,ggparent);
/**/
         intco[i].dihedral=0.0;
         intco[i].connect_type=LC_CONNECTION;
         }
      else
         {
         /* DETERMINE THE SIGN OF THE ANGLE.  THE CROSS OF CROSS1 AND CROSS2 */
         /* LIES IN THE SAME OR OPPOSITE DIRECTION OF BOND2.  IF SAME, ANGLE */
         /* IS POSITIVE, ELSE ANGLE IS NEGATIVE. */
         intco[i].dihedral=vec_angle(cross1,cross2);
         vec_crossprod(cross1,cross2,cross3);
         if (vec_dotprod(bond2,cross3) < 0.0) intco[i].dihedral *= -1.0;
         intco[i].connect_type=TCT_CONNECTION;
         }
      }

   return(TRUE);
   }