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

#define MAX_MM2_ATOM_TYPES 109

int write_sybyl_mol(outfile)
/*============================================================================*/
/* WRITTEN: J.R.HURST
/* LAST MODIFICATION: 5 FEBRUARY 1993 MVG (ADDED ADDITIONAL ATOM TYPES 61-109)
*/
FILE *outfile;
{
  int bond_cnt;
  int i,j;		/* loop variables */
  int status;		/* RETURN STATUS FROM USER I/O ROUTINES. */
  char mol_name[30];	/* MOLECULE NAME. */

  /* FORM MOLECULE NAME FROM THE INPUT FILE NAME. */
  strcpy(mol_name,input_file_name);
  *strrchr(mol_name,'.')=0;		/* TRIM OFF EXTENSION. */
  mol_name[20] = 0;
  fprintf(outfile,"%4d MOL          %-20s              0\n",atom_cnt,mol_name);

  /* WRITE OUT EACH ATOM AND COUNT THE NUMBER OF BONDS. */
  for(bond_cnt = i = 0 ; i < atom_cnt ; i++)
    {
    fprintf(outfile,"%4d%4d%9.4f%9.4f%9.4f     \n",i+1,sybyl_atom_type(i),
       intco[i].vector[0],intco[i].vector[1],intco[i].vector[2]);

    bond_cnt += intco[i].bond_cnt;
    }

  /* WRITE OUT THE BONDS. */
  fprintf(outfile,"%4d MOL\n",bond_cnt/2);
  for(bond_cnt = 1, i = 0 ; i < atom_cnt ; i++)
    {
    for(j = 0 ; j < intco[i].bond_cnt ; j++,bond_cnt++)
      {
      if(i < intco[i].bonds[j])
        fprintf(outfile,"%4d%4d%4d         %4d\n",bond_cnt,i+1,intco[i].bonds[j]+1,1);
      }
    } 

  /* NO FEATURES. */
  fprintf(outfile,"   0 MOL\n");
  
  return(1);
  }

int sybyl_atom_type(index)
/*============================================================================*/
/* PURPOSE: DETERMINE THE SYBYL ATOM TYPE FOR THE SPECIFIED ATOM.  DEPENDING ON
/*	THE GEOMETRY SOURCE, WE CAN EITHER DO A GOOD JOB (MM2/MM3) OR A CRUDE
/*	JOB (ALL ELSE).
*/
   int index;
   {
   /* NOTE THAT THE FOLLOWING TRANSLATION TABLE IS SUPERCEDED IN SOME SPECIAL */
   /* CASES BY CODE SEGMENTS BELOW; IF NOT A SPECIAL CASE, THE TABLE IS USED. */
   static int sybyl_type[MAX_MM2_ATOM_TYPES] =
      {
      1,2,2,4,13,8,9,5,28,7,
      16,15,14,17,10,10,29,30,27,20,
      13,1,13,13,12,986,986,13,1,2,
      915,931,948,917,933,13,11,2,31,19,
      9,18,6,13,7,19,8,13,8,2,
      984,987,990,918,934,1,2,2,988,12,
      910,910,911,911,906,906,2,4,9,9,
      2,6,13,26,9,9,9,9,9,9,
      9,9,9,9,9,9,9,9,9,9,
      9,9,9,9,9,9,9,9,9,9,
      9,9,9,10,10,4,19,6,19
      };

   if (geometry_source==TYPE_MM2_INPUT || geometry_source==TYPE_MM3_INPUT)
      {
      /* SPECIAL KNOWLEDGE USED FOR MM2/MM3 TYPE SOURCE GEOMETRY. */
      if (intco[index].source_type==2)
         {
         /* SPECIAL CASE: CHECK AROMATICITY OF SP2 CARBON (MM2 TYPE 2). */
         if (intco[index].source_aromatic) return(3);
         else return(2);
         }
      else if (intco[index].source_type==37)
         {
         /* SPECIAL CASE: CHECK AROMATICITY OF NITROGEN (MM2 TYPE 37). */
         if (intco[index].source_aromatic) return(11);
         else return(6);
         }
      else if (intco[index].source_type==41)
         {
         /* SPECIAL CASE: CHECK AROMATICITY OF FURAN OXYGEN (MM2 TYPE 41). */
         if (intco[index].source_aromatic) return(9);
         else return(8);
         }
      else return(sybyl_type[intco[index].source_type-1]);
      }
   else if (geometry_source==TYPE_SYBYL_MOL || geometry_source==TYPE_SYBYL_MOL2)
      {
      /* ALREADY CORRECT IF SYBYL WAS THE INCOMING TYPE. */
      return(intco[index].source_type);
      }
   else
      {
      /* ALL NON MM2/MM3 TYPE SOURCE GEOMETRIES; DO A CRUDE TRANSLATION. */
      if(strcmp(intco[index].atom_name,"C ") == 0) return(1);
      if(strcmp(intco[index].atom_name,"H ") == 0) return(13);
      if(strcmp(intco[index].atom_name,"N ") == 0) return(5);
      if(strcmp(intco[index].atom_name,"O ") == 0) return(8);
      if(strcmp(intco[index].atom_name,"S ") == 0) return(10);
      if(strcmp(intco[index].atom_name,"P ") == 0) return(12);
      if(strcmp(intco[index].atom_name,"F ") == 0) return(16);
      if(strcmp(intco[index].atom_name,"Cl") == 0) return(15);
      if(strcmp(intco[index].atom_name,"Br") == 0) return(14);
      if(strcmp(intco[index].atom_name,"I ") == 0) return(17);
      if(strcmp(intco[index].atom_name,"Si") == 0) return(27);
      if(strcmp(intco[index].atom_name,"Na") == 0) return(21);
      if(strcmp(intco[index].atom_name,"Li") == 0) return(24);
      if(strcmp(intco[index].atom_name,"Al") == 0) return(25);
      if(strcmp(intco[index].atom_name,"LP") == 0) return(20);
      if(strcmp(intco[index].atom_name,"Du") == 0) return(26);

      printf("Atom %s not recognized - assumed to be Carbon.\n",intco[index].atom_name);
      return(1);
      }
}

int write_sybyl_mol2(outfile)
/*============================================================================*/
/* PURPOSE: WRITE A SYBYL .MOL2 FILE.  NOTE WE ONLY HANDLE THREE TYPES OF
/*	RECORDS IN THE .MOL2 FILE (MOLECULE, ATOM, AND BOND).
*/
   FILE *outfile;
   {
   int bond_cnt;	/* NUMBER OF BONDS IN MOLECULE. */
   int natom;		/* LOOP INDEX; WHICH ATOM. */
   int nbond;		/* LOOP INDEX; WHICH BOND. */
   int natom_bond;	/* LOOP INDEX; WHICH BOND FOR NATOM. */

   char *sybyl_mnemonic_type();

   /* MAKE A QUICK COUNT OF THE NUMBER OF BONDS IN THE MOLECULE. */
   bond_cnt=0;
   for (natom=0; natom<atom_cnt; natom++)
      {
      bond_cnt+=intco[natom].bond_cnt;
      }
   bond_cnt/=2;

   /* CREATE THE MOLECULE SECTION. */
   fprintf(outfile,"@<TRIPOS>MOLECULE\n");
   fprintf(outfile,"%s\n",source_title);		/* NAME. */
   fprintf(outfile,"%d %d 0 0 0\n",atom_cnt,bond_cnt);	/* NATOMS, NBONDS. */
   fprintf(outfile,"SMALL\n");				/* MOL TYPE. */
   fprintf(outfile,"NO_CHARGES\n");			/* CHARGE TYPE. */
   fprintf(outfile,"****\n");				/* STATUS BITS. */
   fprintf(outfile,"%s\n",header_text[0]);		/* COMMENT. */

   /* CREATE THE ATOM SECTION. */
   /* NOTE THAT THE MNEMONIC IS ONLY CORRECT IF THE SOURCE WAS SYBYL!! */
   fprintf(outfile,"@<TRIPOS>ATOM\n");
   for (natom=0; natom<atom_cnt; natom++)
      {
      fprintf(outfile,"%d %s %f %f %f %s\n",natom+1,intco[natom].atom_name,
         intco[natom].vector[0],intco[natom].vector[1],intco[natom].vector[2],
         sybyl_mnemonic_type(sybyl_atom_type(natom)));
      }

   /* CREATE THE BOND SECTION. */
   fprintf(outfile,"@<TRIPOS>BOND\n");
   nbond=1;
   for (natom=0; natom<atom_cnt; natom++)
      {
      for (natom_bond=0; natom_bond<intco[natom].bond_cnt; natom_bond++)
         {
         if (natom<intco[natom].bonds[natom_bond])
            fprintf(outfile,"%d %d %d 1\n",nbond++,natom+1,intco[natom].bonds[natom_bond]+1);
         }
      }

   return(TRUE);
   }