/*============================================================================*/
/* FILENAME: MM2_PARAMS.C (MM2_PARAMS)
/* USAGE: MM2_PARAMS [-A] [-D] {MM2/MM3 OUTPUT FILE WITH ERRORS}
/*	-A	AUTOMATICALLY TAKE ALL DERIVED PARAMETERS.
/*	-D	PRINTS OUT DEBUGGING INFO AS TO HOW PARAMETERS ARE DERIVED.
/* PURPOSE: SIFT THROUGH MM2/MM3 OUTPUT FILE FOR ERRORS REGARDING MISSING
/*	PARAMETERS AND DERIVE THOSE PARAMETERS BASED ON SIMILARITIES
/*	BETWEEN THE INVOLVED ATOM TYPES AND EXISTING PARAMETERS FOR OTHER
/*	ATOM TYPES (OR BY ASKING USER IF NO SIMILAR CASES EXIST).
/*	NOTE THAT THERE ARE SIGNIFICANT DIFFERENCES IN THE DETAILS OF HANDLING
/*	MM2 VERSUS MM3 OUTPUT FILES (DIFFERENT BUILTIN PARAMETER FILE FORMATS,
/*	DIFFERENT SITE PARAMETER FILES, DIFFERENT FILENAMES, DIFFERENT ERROR
/*	MESSAGE FORMATS, ETC).
/* REFERENCE: JOURNAL OF COMPUTATIONAL CHEMISTRY, VOLUME 12, NUMBER 7,
/*	844-849 (1991) DORA SCHNUR, MARK GRIESHABER, J. PHILLIP BOWEN,
/*	"DEVELOPMENT OF AN INTERNAL SEARCHING ALGORITHM FOR PARAMETERIZATION
/*	OF THE MM2/MM3 FORCE FIELDS".
/* COPYRIGHT 1992 MONSANTO COMPANY
/* WRITTEN: M.V.GRIESHABER + D.M.SCHNUR
/* LAST MODIFICATION: 6 APRIL 1992 MVG (ADDED COPYRIGHT NOTICE)
/*============================================================================*/
#include "utility.h"

#define MAX_TORSION_ERRORS 100		/* MAX ALLOWED TORSION ERRORS. */
#define MAX_BOND_ERRORS 100		/* MAX ALLOWED BOND ERRORS. */
#define MAX_ANGLE_ERRORS 100		/* MAX ALLOWED (2D) ANGLE ERRORS. */
#define MAX_OUTOFPLANE_ERRORS 100	/* MAX ALLOWED OUT-OF-PLANE ERRORS. */

#define MAX_TORSION_PARAMS_NORMAL 900	/* MAX ALLOWED NORMAL TORSION PARAMETERS. */
#define MAX_TORSION_PARAMS_4RING 200	/* MAX ALLOWED 4 MEMBERED RING TORSION PARAMETERS. */
#define MAX_TORSION_PARAMS_ALLENE 20	/* MAX ALLOWED ALLENE TORSION PARAMETERS. */
#define MAX_BOND_PARAMS 400		/* MAX ALLOWED BOND PARAMETERS. */
#define MAX_ANGLE_PARAMS_NORMAL 600	/* MAX ALLOWED NORMAL (2D) ANGLE PARAMETERS. */
#define MAX_ANGLE_PARAMS_4RING 100	/* MAX ALLOWED 4 MEMBERED RING (2D) ANGLE PARAMETERS. */
#define MAX_ANGLE_PARAMS_3RING 100	/* MAX ALLOWED 3 MEMBERED RING (2D) ANGLE PARAMETERS. */
#define MAX_OUTOFPLANE_PARAMS 200	/* MAX ALLOWED OUT-OF-PLANE PARAMETERS. */

#define MAX_ATOM_TYPES 60		/* MAX MM2/MM3 ATOM TYPES. */
#define MAX_TRANSFORMS 3		/* MAX TRANSFORMATIONS FOR ANY ATOM TYPE. */
#define MAX_ATOMS 256			/* MAX ATOMS IN MOLECULE STRUCTURE. */
#define MAX_BONDS_PER_ATOM 8		/* MAX BONDS TO ONE ATOM. */

#define MM3_END_OF_SECTION '9'		/* END OF MM3 PARAMETER SECTION FLAG. */

#define NORMAL 1			/* NORMAL (DEFAULT) ENVIRONMENT. */
#define ALLENE 2			/* ALLENE ENVIRONMENT. */
#define THREE_RING 3			/* 3-MEMBERED RING ENVIRONMENT. */
#define FOUR_RING 4			/* 4-MEMBERED RING ENVIRONMENT. */

#define NO_MATCH 1			/* NO SIMILAR PARAMETER FOUND. */
#define EXACT_MATCH 2			/* EXACT PARAMETER FOUND. */
#define SIMILAR_MATCH 3			/* SIMILAR PARAMETER FOUND. */

#define READ_ACCESS 4			/* FLAG FOR ACCESS ROUTINE. */
#define ATOM_TYPE 0			/* ARRAY INDEX FOR NEW ATOM TYPE. */
#define COST 1				/* ARRAY INDEX FOR COST. */

struct torsion_data
   {
   int atom_types[4];	/* ATOM TYPES IN TORSION. */
   int atom_numbers[4];	/* ATOM NUMBERS IN TORSION (ERRORS ONLY). */
   double v1;		/* ONE-FOLD TORSIONAL CONSTANT. */
   double v2;		/* TWO-FOLD TORSIONAL CONSTANT. */
   double v3;		/* THREE-FOLD TORSIONAL CONSTANT. */
   int orig_env;	/* ORIGINAL ENVIRONMENT OF TORSION (ERRORS ONLY). */
   int save_param;	/* FLAG; SAVE THIS PARAMETER (ERRORS ONLY). */
   };

struct bond_data
   {
   int atom_types[2];	/* ATOM TYPES IN BOND. */
   double s;		/* STRETCHING CONSTANT. */
   double l0;		/* MINIMUM ENERGY BOND LENGTH. */
   double slps;		/* SLOPE OF BOND-ORDER STRETCHING EQUATION. */
   double slpt;		/* SLOPE OF BOND-ORDER LENGTH EQUATION. */
   double bmom;		/* BOND MOMENT. */
   int save_param;	/* FLAG; SAVE THIS PARAMETER (ERRORS ONLY). */
   };

struct angle_data
   {
   int atom_types[3];	/* ATOM TYPES IN ANGLE. */
   int atom_numbers[3];	/* ATOM NUMBERS IN ANGLE (ERRORS ONLY). */
   double b;		/* BENDING CONSTANT. */
   double t;		/* MINIMUM ENERGY BOND ANGLE (ALL ENVIRONMENTS). */
   int orig_env;	/* ORIGINAL ENVIRONMENT OF ANGLE (ERRORS ONLY). */
   int save_param;	/* FLAG; SAVE THIS PARAMETER (ERRORS ONLY). */
   };

struct outofplane_data
   {
   int atom_types[3];	/* ATOM TYPES IN OUT OF PLANE. */
   double b;		/* BENDING CONSTANT. */
   int save_param;	/* FLAG; SAVE THIS PARAMETER (ERRORS ONLY). */
   };

/* REMEMBER EACH ERROR FOUND IN THE MM2 OUTPUT FILE. */
struct torsion_data torsion_errors[MAX_TORSION_ERRORS];
struct bond_data bond_errors[MAX_BOND_ERRORS];
struct angle_data angle_errors[MAX_ANGLE_ERRORS];
struct outofplane_data outofplane_errors[MAX_OUTOFPLANE_ERRORS];

int ntorsion_errors;	/* NUMBER OF TORSION (DIHEDRAL) ERRORS. */
int nbond_errors;	/* NUMBER OF BOND ERRORS. */
int nangle_errors;	/* NUMBER OF (2D) ANGLE ERRORS. */
int noutofplane_errors;	/* NUMBER OF OUT-OF-PLANE ERRORS. */

/* REMEMBER ALL KNOWN PARAMETERS FOR COMPARISONS (TORSION AND ANGLE DATA IS */
/* SPECIFIC FOR THE STRUCTURAL ENVIRONMENT). */
/* NOTE: USAGE OF TORSION_PARAMS_ALLENE IS *NOT* CURRENTLY IMPLEMENTED SINCE */
/* THERE IS A BUG IN MM2 SUCH THAT THE ALLENE PARAMETERS DO NOT WORK; HOWEVER */
/* WE LEAVE THE DATA STRUCTURE HERE FOR FUTURE ACTIVATION. */
struct torsion_data torsion_params_normal[MAX_TORSION_PARAMS_NORMAL];
struct torsion_data torsion_params_4ring[MAX_TORSION_PARAMS_4RING];
struct torsion_data torsion_params_allene[MAX_TORSION_PARAMS_ALLENE];
struct bond_data bond_params[MAX_BOND_PARAMS];
struct angle_data angle_params_normal[MAX_ANGLE_PARAMS_NORMAL];
struct angle_data angle_params_4ring[MAX_ANGLE_PARAMS_4RING];
struct angle_data angle_params_3ring[MAX_ANGLE_PARAMS_3RING];
struct outofplane_data outofplane_params[MAX_OUTOFPLANE_PARAMS];

int ntorsion_params_normal;	/* NUMBER OF NORMAL TORSION PARAMETERS. */
int ntorsion_params_4ring;	/* NUMBER OF 4 MEMBER RING TORSION PARAMETERS. */
int ntorsion_params_allene;	/* NUMBER OF ALLENE TORSION PARAMETERS. */
int nbond_params;		/* NUMBER OF BOND PARAMETERS. */
int nangle_params_normal;	/* NUMBER OF NORMAL (2D) ANGLE PARAMETERS. */
int nangle_params_4ring;	/* NUMBER OF 4 MEMBERED RING (2D) ANGLE PARAMETERS. */
int nangle_params_3ring;	/* NUMBER OF 3 MEMBERED RING (2D) ANGLE PARAMETERS. */
int noutofplane_params;		/* NUMBER OF OUT-OF-PLANE PARAMETERS. */
int debug_on;			/* FLAG; EXTRA DEBUGGING INFORMATION PRINTED. */
int auto_take;			/* FLAG; AUTOMATICALLY TAKE ALL DEFAULTS. */

struct				/* MOLECULE STRUCTURE, CONNECTIVITY ONLY. */
	{
	int nbonds;			/* NUMBER OF BONDS TO THIS ATOM. */
	int bonds[MAX_BONDS_PER_ATOM];	/* ATOMS BONDED TO THIS ATOM. */
	} connect[MAX_ATOMS];
int natoms;			/* NUMBER OF ATOMS IN MOLECULE STRUCTURE. */

int main(argc,argv)
/*============================================================================*/
/* PURPOSE: CREATE MISSING MM2/MM3 PARAMETERS BASED ON SIMILARITY TO OTHER
/*	(KNOWN) PARAMETERS.
*/
   int argc;
   char* argv[];
   {
   char program[10];		/* NAME OF PROGRAM USED (MM2 OR MM3). */
   char error_filename[100];	/* NAME OF MM2 OUTPUT FILE CONTAINING ERRORS. */
   char param_filename[100];	/* NAME OF CREATED NEW PARAMETER FILE. */

   int initialize();		/* VERIFY INPUTS, READ KNOWN PARAMETERS. */
   void digest_mm2_errors();	/* COLLECT ERRORS FROM MM2 .OUT FILE. */
   void digest_mm3_errors();	/* COLLECT ERRORS FROM MM3 .OUT FILE. */
   int make_parameters();	/* GENERATE MISSING PARAMETERS. */
   void write_parameters();	/* CREATE THE MM2 ADDITIONAL PARAMETER FILE. */

   /* VALIDATE INPUT AND READ KNOWN PARAMETERS FOR LATER COMPARISON. */
   if (!initialize(argc,argv,program,error_filename,param_filename)) return(PROGRAM_FAILED);

   /* PICK THE ERRORS OUT OF THE MM2/MM3 OUTPUT FILE. */
   if (strcmp(program,"mm2")==0)
      digest_mm2_errors(error_filename);
   else
      digest_mm3_errors(error_filename);

   /* CREATE THE MISSING PARAMETERS BY SIMILARITY OR DIRECT USER QUERY AND */
   /* WRITE OUT THE CORRECTLY FORMATTED MM2 ADDITIONAL INPUT PARAMETER FILE. */
   if (make_parameters())
      {
      write_parameters(param_filename);
      }

   return(PROGRAM_SUCCEEDED);
   }

int initialize(argc,argv,program,error_filename,param_filename)
/*============================================================================*/
/* PURPOSE: VERIFY INPUTS AND FILE EXISTENCE, AND READ KNOWN PARAMETERS
/*	FOR LATER COMPARISONS.
*/
   int argc;
   char* argv[];
   char program[];		/* NAME OF PROGRAM USED (MM2 OR MM3). */
   char error_filename[];	/* NAME OF MM2 OUTPUT FILE CONTAINING ERRORS. */
   char param_filename[];	/* NAME OF CREATED NEW PARAMETER FILE. */
   {
   char struct_filename[100];	/* NAME OF MOLECULE STRUCTURE FILE. */
   char builtin_filename[100];	/* NAME OF BUILTIN DEFAULT PARAMETER FILE. */
   char site_filename[100];	/* NAME OF LOCAL SITE PARAMETER FILE. */
   int option;			/* OPTION CHARACTER. */

   extern int optind;		/* ARGV INDEX OF NEXT OPTION TO PROCESS. */
   extern int opterr;		/* ERROR CONTROL FOR OPTION PROCESSING. */
   extern char *optarg;		/* OPTION ARGUMENT. */
   int getopt();		/* GET OPTION LETTER FROM ARGUMENT VECTOR. */
   char* getenv();		/* GET ENVIRONMENT NAME VALUE. */

   /* CHECK PROPER USAGE. */
   if (argc>4 || argc<2)
      {
      printf("Usage: mm2_params [-a] [-d] {mm2/mm3 output file containing errors}\n");
      printf("       -a = automatically take all parameters found.\n");
      printf("       -d = debug printing on.\n");
      return(FALSE);
      }

   opterr=0;	/* DISABLE ERROR MESSAGE ON BAD OPTION (SIMPLY IGNORE). */
   debug_on=auto_take=FALSE;
   while ((option=getopt(argc,argv,"ad")) != -1)
      {
      if (option=='a') auto_take=TRUE;
      else if (option=='d') debug_on=TRUE;
      }

   /* FORM THE NAME OF THE MM2/MM3 OUTPUT FILE (WHICH CONTAINS THE ERRORS), */
   /* AND DETERMINE THE NAME OF THE ORIGINATING PROGRAM. */
   strcpy(error_filename,argv[optind]);
   new_extension(error_filename,".out");
   if (strstr(error_filename,"_mm2")!=NULL) strcpy(program,"mm2");
   else if (strstr(error_filename,"_mm3")!=NULL) strcpy(program,"mm3");
   else
      {
      printf("Error: Illegal MM2/MM3 output file name \"%s\".\n",error_filename);
      return(FALSE);
      }
   if (access(error_filename,READ_ACCESS)!=0)
      {
      printf("Error: Cannot read MM2/MM3 output file \"%s\".\n",error_filename);
      return(FALSE);
      }

   /* GET THE MOLECULE STRUCTURE FROM THE ORIGINAL MM2/MM3 INPUT FILE. */
   strcpy(struct_filename,error_filename);
   *strrchr(struct_filename,'_')=0;	/* TRUNCATE "_MM2.OUT". */
   if (strcmp(program,"mm2")==0)
      new_extension(struct_filename,".mm2");
   else
      new_extension(struct_filename,".mm3");
   if (!read_mm2_input(struct_filename)) return(FALSE);

   /* READ THE BUILTIN DEFAULT (CAME WITH MM2 PROGRAM) PARAMETERS. */
   strcpy(builtin_filename,getenv("LOCAL_DATA"));
   if (strcmp(program,"mm2")==0)
      {
      strcat(builtin_filename,"/mm2.builtin_parameters");
      if (!read_default_mm2_params(builtin_filename)) return(FALSE);
      }
   else
      {
      strcat(builtin_filename,"/KONST89.MM3");
      if (!read_default_mm3_params(builtin_filename)) return(FALSE);
      }

   /* READ THE LOCAL SITE (LOCALLY OPTIMIZED) PARAMETERS. */
   strcpy(site_filename,getenv("LOCAL_DATA"));
   if (strcmp(program,"mm2")==0)
      strcat(site_filename,"/mm2.site_parameters");
   else
      strcat(site_filename,"/mm3.site_parameters");
   if (!read_local_params(site_filename)) return(FALSE);

   /* FORM THE NAME OF THE NEW INPUT PARAMETER FILE TO BE CREATED. */
   strcpy(param_filename,error_filename);
   *strrchr(param_filename,'_')=0;	/* TRUNCATE "_MM2.OUT". */
   new_extension(param_filename,".para");

   return(TRUE);
   }

int read_mm2_input(struct_filename)
/*============================================================================*/
/* PURPOSE: READ IN MOLECULE STRUCTURE (CONNECTIVITY) FROM AN MM2 INPUT FILE.
*/
   char struct_filename[];	/* NAME OF MM2 INPUT (STRUCTURE) FILE. */
   {
   FILE* struct_file;		/* OPENED MM2 INPUT FILE. */
   char line[130];		/* RAW LINE FROM MM2 INPUT FILE. */
   int method;			/* FLAG; WHETHER PI SYSTEM INFO PRESENT. */
   int nconnected_atom_lists;	/* NUMBER OF CONNECTED ATOM LISTS. */
   int nattached_atoms;		/* NUMBER OF ATTACHED ATOMS. */
   int i;			/* LOOP INDEX. */

   /* OPEN THE SPECIFIED MM2 INPUT STRUCTURE FILE. */
   if ((struct_file=fopen(struct_filename,"r"))==NULL)
      {
      printf("Error: Cannot read MM2/MM3 input file \"%s\".\n",struct_filename);
      return(FALSE);
      }

   /* GET THE METHOD AND NUMBER OF ATOMS. */
   if (fgets(line,sizeof line,struct_file)==NULL)
      return(errmsg("Reading MM2/MM3 input file (number of atoms)",FALSE));
   sscanf(&line[60],"%1d%4d",&method,&natoms);

   /* SKIP THE AROMATIC FLAGS (THIS MAY CONTINUE ACROSS MULTIPLE CARDS). */
   if (method)
      {
      do
         {
         if (fgets(line,sizeof line,struct_file)==NULL)
            return(errmsg("Reading MM2/MM3 input file (aromatic flags)",FALSE));
         } while (line[79]=='1');
      }

   /* GET NUMBER OF CONNECTED ATOM LISTS AND NUMBER OF ATTACHED ATOMS. */
   if (fgets(line,sizeof line,struct_file)==NULL)
      return(errmsg("Reading MM2/MM3 input file (number of atom lists)",FALSE));
   sscanf(line,"%5d %*f %5d",&nconnected_atom_lists,&nattached_atoms);

   /* BUILD THE BONDS BY READING IN THE CONNECTED AND ATTACHED ATOM LISTS. */
   for (i=0; i<nconnected_atom_lists; i++)
      {
      if (fgets(line,sizeof line,struct_file)==NULL)
         return(errmsg("Reading MM2/MM3 input file (connected atom lists)",FALSE));
      if (!process_connected_atom_list(line)) return(FALSE);
      }
   for (i=0; i<nattached_atoms; i+=8)
      {
      if (fgets(line,sizeof line,struct_file)==NULL)
         return(errmsg("Reading MM2/MM3 input file (attached atoms)",FALSE));
      if (!process_attached_atom_list(line)) return(FALSE);
      }

   fclose(struct_file);
   return(TRUE);
   }

int process_connected_atom_list(line)
/*============================================================================*/
/* PURPOSE: PROCESS AN MM2 CONNECTED ATOM LIST INTO THE INDIVIDUAL BONDS.
/*	NOTE THAT THE SECOND ATOM OF THE PREVIOUS PAIR IS THE FIRST ATOM OF
/*	THE NEXT PAIR.
*/
   char line[];		/* RAW CONNECTED ATOM LIST FROM MM2 INPUT FILE. */
   {
   int i;		/* INDEX OF FIRST ATOM OF PAIR IN LINE. */
   int first_atom;	/* FIRST ATOM OF VALID PAIR. */
   int second_atom;	/* SECOND ATOM OF VALID PAIR. */

   /* TAKE APART A CONNECTED ATOM LIST AN ATOM PAIR AT A TIME. */
   i=0;
   while (i<15 && atoi(&line[(i+1)*5])!=0)
      {
      /* GET THE VALID ATOM PAIR. */
      sscanf(&line[i*5],"%5d%5d",&first_atom,&second_atom);
      first_atom--;	/* ADJUST TO INDEX FROM 0 (MM2 COUNTS FROM 1). */
      second_atom--;
      i++;		/* MOVE TO NEXT ATOM PAIR. */

      /* ADD THE BONDS TO EACH ATOM'S RESPECTIVE BOND LIST. */
      if (connect[first_atom].nbonds>=MAX_BONDS_PER_ATOM || connect[second_atom].nbonds>=MAX_BONDS_PER_ATOM)
         return(errmsg("Too many bonds to a single atom",FALSE));
      connect[first_atom].bonds[connect[first_atom].nbonds++] = second_atom;
      connect[second_atom].bonds[connect[second_atom].nbonds++] = first_atom;
      }

   return(TRUE);
   }

int process_attached_atom_list(line)
/*============================================================================*/
/* PURPOSE: PROCESS AN MM2 ATTACHED ATOM LIST INTO THE INDIVIDUAL BONDS.
/*	NOTE THAT EACH PAIR IS COMPLETELY INDEPENDENT.
*/
   char line[];		/* RAW ATTACHED ATOM LIST FROM MM2 INPUT FILE. */
   {
   int i;		/* INDEX OF FIRST ATOM OF PAIR IN LINE. */
   int first_atom;	/* FIRST ATOM OF VALID PAIR. */
   int second_atom;	/* SECOND ATOM OF VALID PAIR. */

   /* TAKE APART AN ATTACHED ATOM LIST AN ATOM PAIR AT A TIME. */
   i=0;
   while (i<15 && atoi(&line[i*5])!=0)
      {
      /* GET THE VALID ATOM PAIR. */
      sscanf(&line[i*5],"%5d%5d",&first_atom,&second_atom);
      first_atom--;	/* ADJUST TO INDEX FROM 0 (MM2 COUNTS FROM 1). */
      second_atom--;
      i+=2;		/* FINISHED THIS PAIR, MOVE AHEAD TO NEXT. */

      /* ADD THE BONDS TO EACH ATOM'S RESPECTIVE BOND LIST. */
      if (connect[first_atom].nbonds>=MAX_BONDS_PER_ATOM || connect[second_atom].nbonds>=MAX_BONDS_PER_ATOM)
         return(errmsg("Too many bonds to a single atom",FALSE));
      connect[first_atom].bonds[connect[first_atom].nbonds++] = second_atom;
      connect[second_atom].bonds[connect[second_atom].nbonds++] = first_atom;
      }

   return(TRUE);
   }

int read_default_mm2_params(builtin_filename)
/*============================================================================*/
/* PURPOSE: READ IN ALL NEEDED DEFAULT PARAMETERS FOR MM2.
*/
   char builtin_filename[];	/* NAME OF THE BUILTIN PARAMETER FILE. */
   {
   FILE* builtin_file;		/* FILE DESCRIPTOR FOR KNOWN PARAMETER FILE. */
   char line[100];		/* RAW LINE FROM PARAMETER FILE. */

   builtin_file=fopen(builtin_filename,"r");

   /* FIND THE START OF THE REGULAR (NON SPECIAL) TORSIONAL PARAMETERS, */
   /* THEN SUCK UP REGULAR TORSION PARAMETERS UNTIL END OF SECTION. */
   find_line(builtin_file,"TORSIONAL PARAMETERS",line,sizeof line);
   find_line(builtin_file,"ANGLE           V1",line,sizeof line);
   while (!find_string(fgets(line,sizeof line,builtin_file),"4-MEMBERED RING"))
      {
      /* SKIP BLANK, FORTRAN CONTROL, AND EXTRA COLUMN HEADER LINES. */
      if (strlen(line)>1 && line[0]!='+' && line[0]!='1' && line[0]!='0' && !find_string(line,"ANGLE"))
         {
         sscanf(line,"%*d %d-%d-%d-%d %lf %lf %lf",
            &torsion_params_normal[ntorsion_params_normal].atom_types[0],
            &torsion_params_normal[ntorsion_params_normal].atom_types[1],
            &torsion_params_normal[ntorsion_params_normal].atom_types[2],
            &torsion_params_normal[ntorsion_params_normal].atom_types[3],
            &torsion_params_normal[ntorsion_params_normal].v1,
            &torsion_params_normal[ntorsion_params_normal].v2,
            &torsion_params_normal[ntorsion_params_normal].v3);
         ntorsion_params_normal++;
         }
      }

   /* NOW AT START OF 4 MEMBERED RING TORSION PARAMETERS, SKIP TO DATA, */
   /* THEN SUCK UP 4 MEMBER TORSIONAL PARAMETERS UNTIL END OF SECTION. */
   find_line(builtin_file,"ANGLE           V1",line,sizeof line);
   while (!find_string(fgets(line,sizeof line,builtin_file),"ALLENES"))
      {
      /* SKIP BLANK, FORTRAN CONTROL, AND EXTRA COLUMN HEADER LINES. */
      if (strlen(line)>1 && line[0]!='+' && line[0]!='1' && line[0]!='0' && !find_string(line,"ANGLE"))
         {
         sscanf(line,"%*d %d-%d-%d-%d %lf %lf %lf",
            &torsion_params_4ring[ntorsion_params_4ring].atom_types[0],
            &torsion_params_4ring[ntorsion_params_4ring].atom_types[1],
            &torsion_params_4ring[ntorsion_params_4ring].atom_types[2],
            &torsion_params_4ring[ntorsion_params_4ring].atom_types[3],
            &torsion_params_4ring[ntorsion_params_4ring].v1,
            &torsion_params_4ring[ntorsion_params_4ring].v2,
            &torsion_params_4ring[ntorsion_params_4ring].v3);
         ntorsion_params_4ring++;
         }
      }

   /* NOW AT START OF ALLENES TORSION PARAMETERS, SKIP TO DATA, THEN */
   /* SUCK UP ALLENES TORSIONAL PARAMETERS UNTIL END OF SECTION. */
   find_line(builtin_file,"ANGLE           V1",line,sizeof line);
   while (!find_string(fgets(line,sizeof line,builtin_file),"MM2 (1987)"))
      {
      /* SKIP BLANK, FORTRAN CONTROL, AND EXTRA COLUMN HEADER LINES. */
      if (strlen(line)>1 && line[0]!='+' && line[0]!='1' && line[0]!='0' && !find_string(line,"ANGLE"))
         {
         sscanf(line,"%*d %d-%d-%d-%d %lf %lf %lf",
            &torsion_params_allene[ntorsion_params_allene].atom_types[0],
            &torsion_params_allene[ntorsion_params_allene].atom_types[1],
            &torsion_params_allene[ntorsion_params_allene].atom_types[2],
            &torsion_params_allene[ntorsion_params_allene].atom_types[3],
            &torsion_params_allene[ntorsion_params_allene].v1,
            &torsion_params_allene[ntorsion_params_allene].v2,
            &torsion_params_allene[ntorsion_params_allene].v3);
         ntorsion_params_allene++;
         }
      }

   /* FIND THE START OF THE BOND PARAMETERS, THEN */
   /* SUCK UP BOND PARAMETERS UNTIL END OF SECTION. */
   find_line(builtin_file,"STRETCHING PARAMETERS",line,sizeof line);
   find_line(builtin_file,"BOND         KS",line,sizeof line);
   while (!find_string(fgets(line,sizeof line,builtin_file),"4-MEMBERED RING"))
      {
      /* SKIP BLANK, FORTRAN CONTROL, AND EXTRA COLUMN HEADER LINES. */
      if (strlen(line)>1 && line[0]!='+' && line[0]!='1' && line[0]!='0' && !find_string(line,"BOND"))
         {
         sscanf(line,"%*d %d-%d %lf %lf %lf %lf %lf",
            &bond_params[nbond_params].atom_types[0],
            &bond_params[nbond_params].atom_types[1],
            &bond_params[nbond_params].s,
            &bond_params[nbond_params].l0,
            &bond_params[nbond_params].bmom,
            &bond_params[nbond_params].slps,
            &bond_params[nbond_params].slpt);

         /* BONDS MARKED WITH 999 ARE PLACE HOLDERS AND ARE INVALID. */
         if (bond_params[nbond_params].s < 999.0) nbond_params++;
         }
      }

   /* FIND THE START OF THE REGULAR (NON SPECIAL) ANGLE PARAMETERS, THEN */
   /* SUCK UP ANGLE PARAMETERS UNTIL END OF SECTION. */
   find_line(builtin_file,"BENDING PARAMETERS",line,sizeof line);
   find_line(builtin_file,"ANGLE         KS",line,sizeof line);
   find_line(builtin_file,"(-CR2-)",line,sizeof line);
   while (!find_string(fgets(line,sizeof line,builtin_file),"TYPE-2 BENDING KS"))
      {
      /* SKIP BLANK, FORTRAN CONTROL, AND EXTRA COLUMN HEADER LINES. */
      if (strlen(line)>1 && line[0]!='+' && line[0]!='1' && line[0]!='0' &&
         !find_string(line,"ANGLE") && !find_string(line,"(-CR2-)"))
         {
         sscanf(line,"%*d %d-%d-%d %lf %lf",
            &angle_params_normal[nangle_params_normal].atom_types[0],
            &angle_params_normal[nangle_params_normal].atom_types[1],
            &angle_params_normal[nangle_params_normal].atom_types[2],
            &angle_params_normal[nangle_params_normal].b,
            &angle_params_normal[nangle_params_normal].t);

         /* ANGLES MARKED WITH 99 ARE PLACE HOLDERS AND ARE INVALID. */
         if (angle_params_normal[nangle_params_normal].b < 99.0) nangle_params_normal++;
         }
      }

   /* FIND THE START OF THE 4-MEMBERED RING ANGLE PARAMETERS, THEN */
   /* SUCK UP 4-MEMBERED RING ANGLE PARAMETERS UNTIL END OF SECTION. */
   find_line(builtin_file,"4-MEMBERED RING",line,sizeof line);
   find_line(builtin_file,"ANGLE         KS",line,sizeof line);
   find_line(builtin_file,"(-CR2-)",line,sizeof line);
   while (!find_string(fgets(line,sizeof line,builtin_file),"TYPE-2 BENDING KS"))
      {
      /* SKIP BLANK, FORTRAN CONTROL, AND EXTRA COLUMN HEADER LINES. */
      if (strlen(line)>1 && line[0]!='+' && line[0]!='1' && line[0]!='0' &&
         !find_string(line,"ANGLE") && !find_string(line,"(-CR2-)"))
         {
         sscanf(line,"%*d %d-%d-%d %lf %lf",
            &angle_params_4ring[nangle_params_4ring].atom_types[0],
            &angle_params_4ring[nangle_params_4ring].atom_types[1],
            &angle_params_4ring[nangle_params_4ring].atom_types[2],
            &angle_params_4ring[nangle_params_4ring].b,
            &angle_params_4ring[nangle_params_4ring].t);
         nangle_params_4ring++;
         }
      }

   /* FIND THE START OF THE 3-MEMBERED RING ANGLE PARAMETERS, THEN */
   /* SUCK UP 3-MEMBERED RING ANGLE PARAMETERS UNTIL END OF SECTION. */
   find_line(builtin_file,"3-MEMBERED RING",line,sizeof line);
   while (!find_string(fgets(line,sizeof line,builtin_file),"TYPE-2 BENDING KS"))
      {
      /* SKIP BLANK, FORTRAN CONTROL, AND EXTRA COLUMN HEADER LINES. */
      if (strlen(line)>1 && line[0]!='+' && line[0]!='1' && line[0]!='0' &&
         !find_string(line,"ANGLE") && !find_string(line,"(-CR2-)"))
         {
         sscanf(line,"%*d %d-%d-%d %lf %lf",
            &angle_params_3ring[nangle_params_3ring].atom_types[0],
            &angle_params_3ring[nangle_params_3ring].atom_types[1],
            &angle_params_3ring[nangle_params_3ring].atom_types[2],
            &angle_params_3ring[nangle_params_3ring].b,
            &angle_params_3ring[nangle_params_3ring].t);
         nangle_params_3ring++;
         }
      }

   /* FIND THE START OF THE OUT OF PLANE PARAMETERS, THEN */
   /* SUCK UP OUT OF PLANE PARAMETERS UNTIL END OF SECTION. */
   find_line(builtin_file,"OUT-OF-PLANE BENDING PARAMETERS",line,sizeof line);
   while (!find_string(fgets(line,sizeof line,builtin_file),"STRETCH-BEND"))
      {
      /* SKIP BLANK AND FORTRAN CONTROL LINES. */
      if (strlen(line)>1 && line[0]!='+' && line[0]!='1' && line[0]!='0')
         {
         outofplane_params[noutofplane_params].atom_types[0]=0;
         sscanf(line,"%*d %d-%d %lf",
            &outofplane_params[noutofplane_params].atom_types[1],
            &outofplane_params[noutofplane_params].atom_types[2],
            &outofplane_params[noutofplane_params].b);
         noutofplane_params++;
         }
      }

   /* DONE READING BUILT IN DEFAULT PARAMETERS, CLOSE THE FILE AND RETURN. */
   fclose(builtin_file);
   return(TRUE);
   }

int read_default_mm3_params(builtin_filename)
/*============================================================================*/
/* PURPOSE: READ IN ALL NEEDED DEFAULT PARAMETERS FOR MM3.  NOTE SLIGHT
/*	TRICKERY WITH ATOM TYPES - THEY ARE PACKED ########, BUT THE FIRST
/*	DIGIT MAY BE BLANK (WHICH PRECLUDES USING SSCANF DIRECTLY).
*/
   char builtin_filename[];	/* NAME OF THE BUILTIN PARAMETER FILE. */
   {
   FILE* builtin_file;		/* FILE DESCRIPTOR FOR KNOWN PARAMETER FILE. */
   char line[100];		/* RAW LINE FROM PARAMETER FILE. */
   char packed_atom_types[20];	/* ATOM TYPES PACKED AS ########. */
   char tmp_str[20];		/* TEMP STRING FOR STRNCPY. */

   builtin_file=fopen(builtin_filename,"r");

   /* FIND THE START OF THE REGULAR (NON SPECIAL) TORSIONAL PARAMETERS, */
   /* THEN SUCK UP REGULAR TORSION PARAMETERS UNTIL END OF SECTION. */
   find_line(builtin_file,"TORSIONAL PARAMETERS",line,sizeof line);
   find_line(builtin_file,"--------------------",line,sizeof line);
   find_line(builtin_file,"--------------------",line,sizeof line);
   while (fgets(line,sizeof line,builtin_file) && line[0]!=MM3_END_OF_SECTION)
      {
      sscanf(line,"%*d %8c %lf %lf %lf",
         packed_atom_types,
         &torsion_params_normal[ntorsion_params_normal].v1,
         &torsion_params_normal[ntorsion_params_normal].v2,
         &torsion_params_normal[ntorsion_params_normal].v3);
      torsion_params_normal[ntorsion_params_normal].atom_types[0]=atoi(strncpy(tmp_str,packed_atom_types,2));
      torsion_params_normal[ntorsion_params_normal].atom_types[1]=atoi(strncpy(tmp_str,&packed_atom_types[2],2));
      torsion_params_normal[ntorsion_params_normal].atom_types[2]=atoi(strncpy(tmp_str,&packed_atom_types[4],2));
      torsion_params_normal[ntorsion_params_normal].atom_types[3]=atoi(strncpy(tmp_str,&packed_atom_types[6],2));
      ntorsion_params_normal++;
      }

   /* SKIP 5 MEMBERED RING SECTION TO GET TO THE START OF THE 4 MEMBERED RING */
   /* TORSION PARAMETERS, SKIP TO DATA, THEN SUCK UP 4 MEMBER TORSIONAL */
   /* PARAMETERS UNTIL END OF SECTION. */
   find_line(builtin_file,"4-MEMBERED RING",line,sizeof line);
   find_line(builtin_file,"--------------------",line,sizeof line);
   while (fgets(line,sizeof line,builtin_file) && line[0]!=MM3_END_OF_SECTION)
      {
      sscanf(line,"%*d %8c %lf %lf %lf",
         packed_atom_types,
         &torsion_params_4ring[ntorsion_params_4ring].v1,
         &torsion_params_4ring[ntorsion_params_4ring].v2,
         &torsion_params_4ring[ntorsion_params_4ring].v3);
      torsion_params_4ring[ntorsion_params_4ring].atom_types[0]=atoi(strncpy(tmp_str,packed_atom_types,2));
      torsion_params_4ring[ntorsion_params_4ring].atom_types[1]=atoi(strncpy(tmp_str,&packed_atom_types[2],2));
      torsion_params_4ring[ntorsion_params_4ring].atom_types[2]=atoi(strncpy(tmp_str,&packed_atom_types[4],2));
      torsion_params_4ring[ntorsion_params_4ring].atom_types[3]=atoi(strncpy(tmp_str,&packed_atom_types[6],2));
      ntorsion_params_4ring++;
      }

   /* FIND START OF ALLENES TORSION PARAMETERS, SKIP TO DATA, THEN */
   /* SUCK UP ALLENES TORSIONAL PARAMETERS UNTIL END OF SECTION. */
   find_line(builtin_file,"ALLENES",line,sizeof line);
   find_line(builtin_file,"--------------------",line,sizeof line);
   while (fgets(line,sizeof line,builtin_file) && line[0]!=MM3_END_OF_SECTION)
      {
      sscanf(line,"%*d %8c %lf %lf %lf",
         packed_atom_types,
         &torsion_params_allene[ntorsion_params_allene].v1,
         &torsion_params_allene[ntorsion_params_allene].v2,
         &torsion_params_allene[ntorsion_params_allene].v3);
      torsion_params_allene[ntorsion_params_allene].atom_types[0]=atoi(strncpy(tmp_str,packed_atom_types,2));
      torsion_params_allene[ntorsion_params_allene].atom_types[1]=atoi(strncpy(tmp_str,&packed_atom_types[2],2));
      torsion_params_allene[ntorsion_params_allene].atom_types[2]=atoi(strncpy(tmp_str,&packed_atom_types[4],2));
      torsion_params_allene[ntorsion_params_allene].atom_types[3]=atoi(strncpy(tmp_str,&packed_atom_types[6],2));
      ntorsion_params_allene++;
      }

   /* FIND THE START OF THE BOND PARAMETERS, THEN */
   /* SUCK UP BOND PARAMETERS UNTIL END OF SECTION. */
   find_line(builtin_file,"BOND PARAMETERS",line,sizeof line);
   find_line(builtin_file,"--------------------",line,sizeof line);
   find_line(builtin_file,"--------------------",line,sizeof line);
   while (fgets(line,sizeof line,builtin_file) && line[0]!=MM3_END_OF_SECTION)
      {
      sscanf(line,"%*d %4c %lf %lf %lf %lf %lf",
         packed_atom_types,
         &bond_params[nbond_params].s,
         &bond_params[nbond_params].l0,
         &bond_params[nbond_params].bmom,
         &bond_params[nbond_params].slps,
         &bond_params[nbond_params].slpt);
      bond_params[nbond_params].atom_types[0]=atoi(strncpy(tmp_str,packed_atom_types,2));
      bond_params[nbond_params].atom_types[1]=atoi(strncpy(tmp_str,&packed_atom_types[2],2));

      /* BONDS MARKED WITH 999 ARE PLACE HOLDERS AND ARE INVALID. */
      if (bond_params[nbond_params].s < 999.0) nbond_params++;
      }

   /* FIND THE START OF THE REGULAR (NON SPECIAL) ANGLE PARAMETERS, THEN */
   /* SUCK UP ANGLE PARAMETERS UNTIL END OF SECTION. */
   find_line(builtin_file,"BENDING PARAMETERS",line,sizeof line);
   find_line(builtin_file,"--------------------",line,sizeof line);
   find_line(builtin_file,"--------------------",line,sizeof line);
   while (fgets(line,sizeof line,builtin_file) && line[0]!=MM3_END_OF_SECTION)
      {
      sscanf(line,"%*d %6c %lf %lf",
         packed_atom_types,
         &angle_params_normal[nangle_params_normal].b,
         &angle_params_normal[nangle_params_normal].t);
      angle_params_normal[nangle_params_normal].atom_types[0]=atoi(strncpy(tmp_str,packed_atom_types,2));
      angle_params_normal[nangle_params_normal].atom_types[1]=atoi(strncpy(tmp_str,&packed_atom_types[2],2));
      angle_params_normal[nangle_params_normal].atom_types[2]=atoi(strncpy(tmp_str,&packed_atom_types[4],2));

      /* ANGLES MARKED WITH 99 ARE PLACE HOLDERS AND ARE INVALID. */
      if (angle_params_normal[nangle_params_normal].b < 99.0) nangle_params_normal++;
      }

   /* FIND THE START OF THE 4-MEMBERED RING ANGLE PARAMETERS, THEN */
   /* SUCK UP 4-MEMBERED RING ANGLE PARAMETERS UNTIL END OF SECTION. */
   find_line(builtin_file,"4-MEMBERED RING",line,sizeof line);
   find_line(builtin_file,"--------------------",line,sizeof line);
   while (fgets(line,sizeof line,builtin_file) && line[0]!=MM3_END_OF_SECTION)
      {
      sscanf(line,"%*d %6c %lf %lf",
         packed_atom_types,
         &angle_params_4ring[nangle_params_4ring].b,
         &angle_params_4ring[nangle_params_4ring].t);
      angle_params_4ring[nangle_params_4ring].atom_types[0]=atoi(strncpy(tmp_str,packed_atom_types,2));
      angle_params_4ring[nangle_params_4ring].atom_types[1]=atoi(strncpy(tmp_str,&packed_atom_types[2],2));
      angle_params_4ring[nangle_params_4ring].atom_types[2]=atoi(strncpy(tmp_str,&packed_atom_types[4],2));
      nangle_params_4ring++;
      }

   /* FIND THE START OF THE 3-MEMBERED RING ANGLE PARAMETERS, THEN */
   /* SUCK UP 3-MEMBERED RING ANGLE PARAMETERS UNTIL END OF SECTION. */
   find_line(builtin_file,"3-MEMBERED RING",line,sizeof line);
   find_line(builtin_file,"--------------------",line,sizeof line);
   while (fgets(line,sizeof line,builtin_file) && line[0]!=MM3_END_OF_SECTION)
      {
      sscanf(line,"%*d %6c %lf %lf",
         packed_atom_types,
         &angle_params_3ring[nangle_params_3ring].b,
         &angle_params_3ring[nangle_params_3ring].t);
      angle_params_3ring[nangle_params_3ring].atom_types[0]=atoi(strncpy(tmp_str,packed_atom_types,2));
      angle_params_3ring[nangle_params_3ring].atom_types[1]=atoi(strncpy(tmp_str,&packed_atom_types[2],2));
      angle_params_3ring[nangle_params_3ring].atom_types[2]=atoi(strncpy(tmp_str,&packed_atom_types[4],2));
      nangle_params_3ring++;
      }

   /* FIND THE START OF THE OUT OF PLANE PARAMETERS, THEN */
   /* SUCK UP OUT OF PLANE PARAMETERS UNTIL END OF SECTION. */
   find_line(builtin_file,"OUT OF PLANE BENDING PARAMETERS",line,sizeof line);
   find_line(builtin_file,"--------------------",line,sizeof line);
   find_line(builtin_file,"--------------------",line,sizeof line);
   while (fgets(line,sizeof line,builtin_file) && line[0]!=MM3_END_OF_SECTION)
      {
      sscanf(line,"%*d %4c %lf",
         packed_atom_types,
         &outofplane_params[noutofplane_params].b);
      outofplane_params[noutofplane_params].atom_types[0]=0;
      outofplane_params[noutofplane_params].atom_types[1]=atoi(strncpy(tmp_str,packed_atom_types,2));
      outofplane_params[noutofplane_params].atom_types[2]=atoi(strncpy(tmp_str,&packed_atom_types[2],2));
      noutofplane_params++;
      }

   /* DONE READING BUILT IN DEFAULT PARAMETERS, CLOSE THE FILE AND RETURN. */
   fclose(builtin_file);
   return(TRUE);
   }

int read_local_params(param_filename)
/*============================================================================*/
/* PURPOSE: READ IN LOCAL SITE PARAMETERS.  THESE ARE PARAMETERS
/*	THAT HAVE BEEN OPTIMIZED AT LOCAL SITE.
*/
   char param_filename[];	/* NAME OF THE LOCAL PARAMETER FILE. */
   {
   FILE* param_file;		/* FILE DESCRIPTOR FOR PARAMETER FILE. */
   char line[100];		/* RAW LINE FROM PARAMETER FILE. */
   int nparameters;		/* NUMBER OF PARAMETERS IN EACH SECTION. */
   int i;			/* LOOP INDEX. */

   /* PERFECTLY VALID NOT TO HAVE ANY LOCAL PARAMETER FILE. */
   if ((param_file=fopen(param_filename,"r"))==NULL) return(TRUE);

   /* READ NORMAL TORSION PARAMETERS. */
   get_data_line(param_file,line,sizeof line);
   sscanf(line,"%d",&nparameters);
   for (i=0; i<nparameters; i++)
      {
      get_data_line(param_file,line,sizeof line);
      sscanf(line,"%d-%d-%d-%d %lf %lf %lf",
         &torsion_params_normal[ntorsion_params_normal].atom_types[0],
         &torsion_params_normal[ntorsion_params_normal].atom_types[1],
         &torsion_params_normal[ntorsion_params_normal].atom_types[2],
         &torsion_params_normal[ntorsion_params_normal].atom_types[3],
         &torsion_params_normal[ntorsion_params_normal].v1,
         &torsion_params_normal[ntorsion_params_normal].v2,
         &torsion_params_normal[ntorsion_params_normal].v3);
      ntorsion_params_normal++;
      }

   /* READ 4 MEMBERED RING TORSION PARAMETERS. */
   get_data_line(param_file,line,sizeof line);
   sscanf(line,"%d",&nparameters);
   for (i=0; i<nparameters; i++)
      {
      get_data_line(param_file,line,sizeof line);
      sscanf(line,"%d-%d-%d-%d %lf %lf %lf",
         &torsion_params_4ring[ntorsion_params_4ring].atom_types[0],
         &torsion_params_4ring[ntorsion_params_4ring].atom_types[1],
         &torsion_params_4ring[ntorsion_params_4ring].atom_types[2],
         &torsion_params_4ring[ntorsion_params_4ring].atom_types[3],
         &torsion_params_4ring[ntorsion_params_4ring].v1,
         &torsion_params_4ring[ntorsion_params_4ring].v2,
         &torsion_params_4ring[ntorsion_params_4ring].v3);
      ntorsion_params_4ring++;
      }

   /* READ BOND PARAMETERS. */
   get_data_line(param_file,line,sizeof line);
   sscanf(line,"%d",&nparameters);
   for (i=0; i<nparameters; i++)
      {
      get_data_line(param_file,line,sizeof line);
      sscanf(line,"%d-%d %lf %lf %lf %lf %lf",
         &bond_params[nbond_params].atom_types[0],
         &bond_params[nbond_params].atom_types[1],
         &bond_params[nbond_params].s,
         &bond_params[nbond_params].l0,
         &bond_params[nbond_params].bmom,
         &bond_params[nbond_params].slps,
         &bond_params[nbond_params].slpt);
      nbond_params++;
      }

   /* READ NORMAL ANGLE PARAMETERS. */
   get_data_line(param_file,line,sizeof line);
   sscanf(line,"%d",&nparameters);
   for (i=0; i<nparameters; i++)
      {
      get_data_line(param_file,line,sizeof line);
      sscanf(line,"%d-%d-%d %lf %lf",
         &angle_params_normal[nangle_params_normal].atom_types[0],
         &angle_params_normal[nangle_params_normal].atom_types[1],
         &angle_params_normal[nangle_params_normal].atom_types[2],
         &angle_params_normal[nangle_params_normal].b,
         &angle_params_normal[nangle_params_normal].t);
      nangle_params_normal++;
      }

   /* READ 4 MEMBERED RING ANGLE PARAMETERS. */
   get_data_line(param_file,line,sizeof line);
   sscanf(line,"%d",&nparameters);
   for (i=0; i<nparameters; i++)
      {
      get_data_line(param_file,line,sizeof line);
      sscanf(line,"%d-%d-%d %lf %lf",
         &angle_params_4ring[nangle_params_4ring].atom_types[0],
         &angle_params_4ring[nangle_params_4ring].atom_types[1],
         &angle_params_4ring[nangle_params_4ring].atom_types[2],
         &angle_params_4ring[nangle_params_4ring].b,
         &angle_params_4ring[nangle_params_4ring].t);
      nangle_params_4ring++;
      }

   /* READ 3 MEMBERED RING ANGLE PARAMETERS. */
   get_data_line(param_file,line,sizeof line);
   sscanf(line,"%d",&nparameters);
   for (i=0; i<nparameters; i++)
      {
      get_data_line(param_file,line,sizeof line);
      sscanf(line,"%d-%d-%d %lf %lf",
         &angle_params_3ring[nangle_params_3ring].atom_types[0],
         &angle_params_3ring[nangle_params_3ring].atom_types[1],
         &angle_params_3ring[nangle_params_3ring].atom_types[2],
         &angle_params_3ring[nangle_params_3ring].b,
         &angle_params_3ring[nangle_params_3ring].t);
      nangle_params_3ring++;
      }

   /* READ OUT OF PLANE PARAMETERS. */
   get_data_line(param_file,line,sizeof line);
   sscanf(line,"%d",&nparameters);
   for (i=0; i<nparameters; i++)
      {
      get_data_line(param_file,line,sizeof line);
      sscanf(line,"%d-%d-%d %lf",
         &outofplane_params[noutofplane_params].atom_types[0],
         &outofplane_params[noutofplane_params].atom_types[1],
         &outofplane_params[noutofplane_params].atom_types[2],
         &outofplane_params[noutofplane_params].b);
      noutofplane_params++;
      }

   /* DONE READING LOCAL PARAMETERS, CLOSE THE FILE AND RETURN. */
   fclose(param_file);
   return(TRUE);
   }

int get_data_line(file,line,max_line_length)
/*============================================================================*/
/* PURPOSE: GET A "DATA" LINE FROM THE SPECIFIED FILE.  THE FILE MUST BE AN
/*	ASCII FILE.  A DATA LINE IS ANY LINE THAT IS NONBLANK AND DOES NOT
/*	START WITH AN OCTOTHORPE (#).  INTERVENING LINES ARE SKIPPED.
*/
   FILE* file;		/* OPENED DATA FILE. */
   char line[];		/* RETRIEVED DATA LINE. */
   int max_line_length;	/* MAXIMUM LENGTH OF LINE. */
   {
   do
      {
      if (fgets(line,max_line_length,file)==NULL) return(FALSE);
      } while(strlen(line)<=1 || line[0]=='#');

   return(TRUE);
   }

void digest_mm2_errors(error_filename)
/*============================================================================*/
/* PURPOSE: DIGEST THE MM2 SPECIFIC OUTPUT FILE CONTAINING ERRORS
/*	INTO THE BASIC ERROR TYPES AND ATOM TYPES.  NOTE THAT ERRORS THAT DO
/*	NOT MATCH THE KNOWN ERROR FORMATS ARE NOT HANDLED BY THIS PROGRAM.
/*	NOTE ALSO THAT FOR TORSIONS AND ANGLES ATOM NUMBERS AS WELL AS TYPES
/*	ARE READ AND STORED (ATOM NUMBERS IN MM2 INDEX FROM 1, NOT 0).
*/
   char error_filename[];	/* NAME OF MM2 OUTPUT FILE WITH ERRORS. */
   {
   FILE* error_file;		/* MM2 OUTPUT FILE CONTAINING ERRORS. */
   char line[100];		/* RAW LINE FROM MM2 OUTPUT FILE. */

   void save_if_new_torsion_error();	/* SAVE ONLY NEW TORSION ERROR. */
   void save_if_new_bond_error();	/* SAVE ONLY NEW BOND ERROR. */
   void save_if_new_angle_error();	/* SAVE ONLY NEW ANGLE ERROR. */
   void save_if_new_oop_error();	/* SAVE ONLY NEW OUT OF PLANE ERROR. */

   ntorsion_errors=nbond_errors=nangle_errors=noutofplane_errors=0;
   error_file=fopen(error_filename,"r");

   /* FIND EACH ERROR LINE IN THE MM2 OUTPUT FILE. */
   while (find_line(error_file," * * * ERROR - CONST",line,sizeof line)!=NULL)
      {
      if (find_string(line,"CONSTS FOR ANG")!=NULL)
         {
         /* EXTRACT ATOM NUMBERS AND TYPES FOR TORSION (DIHEDRAL) ERRORS. */
         sscanf(find_string(line,"ANG"),"ANG %d %d %d %d (TYPE %d %d %d %d",
            &torsion_errors[ntorsion_errors].atom_numbers[0],
            &torsion_errors[ntorsion_errors].atom_numbers[1],
            &torsion_errors[ntorsion_errors].atom_numbers[2],
            &torsion_errors[ntorsion_errors].atom_numbers[3],
            &torsion_errors[ntorsion_errors].atom_types[0],
            &torsion_errors[ntorsion_errors].atom_types[1],
            &torsion_errors[ntorsion_errors].atom_types[2],
            &torsion_errors[ntorsion_errors].atom_types[3]);
         torsion_errors[ntorsion_errors].atom_numbers[0]--;	/* ADJUST TO INDEX FROM 0. */
         torsion_errors[ntorsion_errors].atom_numbers[1]--;
         torsion_errors[ntorsion_errors].atom_numbers[2]--;
         torsion_errors[ntorsion_errors].atom_numbers[3]--;
         torsion_errors[ntorsion_errors].orig_env=get_torsion_environment(torsion_errors[ntorsion_errors].atom_numbers);

         /* SAVE THIS ERROR ONLY IF IT IS A *NEW* ERROR. */
         save_if_new_torsion_error();
         }

      else if (find_string(line,"CONSTANTS FOR BOND")!=NULL)
         {
         /* EXTRACT ATOM TYPES FOR BOND ERRORS. */
         sscanf(find_string(line,"TYPE"),"TYPE %d - %d",
            &bond_errors[nbond_errors].atom_types[0],
            &bond_errors[nbond_errors].atom_types[1]);

         /* SAVE THIS ERROR ONLY IF IT IS A *NEW* ERROR. */
         save_if_new_bond_error();
         }

      else if (find_string(line,"CONSTANTS FOR ANGLE")!=NULL)
         {
         /* EXTRACT ATOM NUMBERS AND TYPES FOR ANGLE (2D) ERRORS. */
         sscanf(find_string(line,"ANGLE"),"ANGLE %d %d %d (TYPE %d %d %d",
            &angle_errors[nangle_errors].atom_numbers[0],
            &angle_errors[nangle_errors].atom_numbers[1],
            &angle_errors[nangle_errors].atom_numbers[2],
            &angle_errors[nangle_errors].atom_types[0],
            &angle_errors[nangle_errors].atom_types[1],
            &angle_errors[nangle_errors].atom_types[2]);
         angle_errors[nangle_errors].atom_numbers[0]--;	/* ADJUST TO INDEX FROM 0. */
         angle_errors[nangle_errors].atom_numbers[1]--;
         angle_errors[nangle_errors].atom_numbers[2]--;
         angle_errors[nangle_errors].orig_env=get_angle_environment(angle_errors[nangle_errors].atom_numbers,angle_errors[nangle_errors].atom_types);

         /* SAVE THIS ERROR ONLY IF IT IS A *NEW* ERROR. */
         save_if_new_angle_error();
         }

      else if (find_string(line,"CONSTANTS FOR OUT-OF-PLANE")!=NULL)
         {
         /* EXTRACT ATOM TYPES FOR OUT OF PLANE ERRORS. */
         sscanf(find_string(line,"TYPE"),"TYPE %d %d %d",
            &outofplane_errors[noutofplane_errors].atom_types[0],
            &outofplane_errors[noutofplane_errors].atom_types[1],
            &outofplane_errors[noutofplane_errors].atom_types[2]);

         /* SAVE THIS ERROR ONLY IF IT IS A *NEW* ERROR. */
         save_if_new_oop_error();
         }
      }

   fclose(error_file);
   return;
   }

void digest_mm3_errors(error_filename)
/*============================================================================*/
/* PURPOSE: DIGEST THE MM3 SPECIFIC OUTPUT FILE CONTAINING ERRORS
/*	INTO THE BASIC ERROR TYPES AND ATOM TYPES.  NOTE THAT ERRORS THAT DO
/*	NOT MATCH THE KNOWN ERROR FORMATS ARE NOT HANDLED BY THIS PROGRAM.
/*	NOTE ALSO THAT FOR TORSIONS AND ANGLES ATOM NUMBERS AS WELL AS TYPES
/*	ARE READ AND STORED (ATOM NUMBERS IN MM3 INDEX FROM 1, NOT 0).
*/
   char error_filename[];	/* NAME OF MM3 OUTPUT FILE WITH ERRORS. */
   {
   FILE* error_file;		/* MM3 OUTPUT FILE CONTAINING ERRORS. */
   char line[100];		/* RAW LINE FROM MM3 OUTPUT FILE. */

   void save_if_new_torsion_error();	/* SAVE ONLY NEW TORSION ERROR. */
   void save_if_new_bond_error();	/* SAVE ONLY NEW BOND ERROR. */
   void save_if_new_angle_error();	/* SAVE ONLY NEW ANGLE ERROR. */
   void save_if_new_oop_error();	/* SAVE ONLY NEW OUT OF PLANE ERROR. */

   ntorsion_errors=nbond_errors=nangle_errors=noutofplane_errors=0;
   error_file=fopen(error_filename,"r");

   /* FIND EACH ERROR LINE IN THE MM3 OUTPUT FILE. */
   while (find_line(error_file," MISSING --- CONSTS",line,sizeof line)!=NULL)
      {
      if (find_string(line,"BOND")!=NULL)
         {
         /* EXTRACT ATOM TYPES FOR BOND ERRORS. */
         sscanf(find_string(line,"TYPE"),"TYPE %d- %d",
            &bond_errors[nbond_errors].atom_types[0],
            &bond_errors[nbond_errors].atom_types[1]);

         /* SAVE THIS ERROR ONLY IF IT IS A *NEW* ERROR. */
         save_if_new_bond_error();
         }

      else if (find_string(line,"JH =")!=NULL)
         {
         /* EXTRACT ATOM NUMBERS AND TYPES FOR ANGLE (2D) ERRORS. */
         sscanf(find_string(line,"FOR"),"FOR %d %d %d (TYPE %d- %d- %d",
            &angle_errors[nangle_errors].atom_numbers[0],
            &angle_errors[nangle_errors].atom_numbers[1],
            &angle_errors[nangle_errors].atom_numbers[2],
            &angle_errors[nangle_errors].atom_types[0],
            &angle_errors[nangle_errors].atom_types[1],
            &angle_errors[nangle_errors].atom_types[2]);
         angle_errors[nangle_errors].atom_numbers[0]--;	/* ADJUST TO INDEX FROM 0. */
         angle_errors[nangle_errors].atom_numbers[1]--;
         angle_errors[nangle_errors].atom_numbers[2]--;
         angle_errors[nangle_errors].orig_env=get_angle_environment(angle_errors[nangle_errors].atom_numbers,angle_errors[nangle_errors].atom_types);

         /* SAVE THIS ERROR ONLY IF IT IS A *NEW* ERROR. */
         save_if_new_angle_error();
         }

      else if (find_string(line,"OUT-OF-PLN")!=NULL)
         {
         /* EXTRACT ATOM TYPES FOR OUT OF PLANE ERRORS. */
         outofplane_errors[noutofplane_errors].atom_types[0]=0;
         sscanf(find_string(line,"TYPE"),"TYPE %d- %d",
            &outofplane_errors[noutofplane_errors].atom_types[1],
            &outofplane_errors[noutofplane_errors].atom_types[2]);

         /* SAVE THIS ERROR ONLY IF IT IS A *NEW* ERROR. */
         save_if_new_oop_error();
         }

      else
         {
         /* EXTRACT ATOM NUMBERS AND TYPES FOR TORSION (DIHEDRAL) ERRORS. */
         sscanf(find_string(line,"FOR"),"FOR %d- %d- %d- %d (TYPE %d- %d- %d- %d",
            &torsion_errors[ntorsion_errors].atom_numbers[0],
            &torsion_errors[ntorsion_errors].atom_numbers[1],
            &torsion_errors[ntorsion_errors].atom_numbers[2],
            &torsion_errors[ntorsion_errors].atom_numbers[3],
            &torsion_errors[ntorsion_errors].atom_types[0],
            &torsion_errors[ntorsion_errors].atom_types[1],
            &torsion_errors[ntorsion_errors].atom_types[2],
            &torsion_errors[ntorsion_errors].atom_types[3]);
         torsion_errors[ntorsion_errors].atom_numbers[0]--;	/* ADJUST TO INDEX FROM 0. */
         torsion_errors[ntorsion_errors].atom_numbers[1]--;
         torsion_errors[ntorsion_errors].atom_numbers[2]--;
         torsion_errors[ntorsion_errors].atom_numbers[3]--;
         torsion_errors[ntorsion_errors].orig_env=get_torsion_environment(torsion_errors[ntorsion_errors].atom_numbers);

         /* SAVE THIS ERROR ONLY IF IT IS A *NEW* ERROR. */
         save_if_new_torsion_error();
         }
      }

   fclose(error_file);
   return;
   }

void save_if_new_torsion_error()
/*============================================================================*/
/* PURPOSE: DETERMINE IT THE TORSION ERROR JUST ADDED TO THE END OF THE LIST
/*	IS A NEW ERROR.  IF SO, SAVE THE ERROR BY INCREMENTING THE NUMBER OF
/*	KNOWN TORSION ERRORS.
*/
   {
   int already_exists;	/* FLAG; T IF TORSION ERROR ALREADY IN LIST. */
   int i;		/* LOOP INDEX; WHICH TORSION ERROR. */

   already_exists=FALSE;
   for (i=0; i<ntorsion_errors; i++)
      {
      if ((torsion_errors[ntorsion_errors].atom_types[0]==torsion_errors[i].atom_types[0] &&
         torsion_errors[ntorsion_errors].atom_types[1]==torsion_errors[i].atom_types[1] &&
         torsion_errors[ntorsion_errors].atom_types[2]==torsion_errors[i].atom_types[2] &&
         torsion_errors[ntorsion_errors].atom_types[3]==torsion_errors[i].atom_types[3] &&
         torsion_errors[ntorsion_errors].orig_env==torsion_errors[i].orig_env) ||
         (torsion_errors[ntorsion_errors].atom_types[0]==torsion_errors[i].atom_types[3] &&
         torsion_errors[ntorsion_errors].atom_types[1]==torsion_errors[i].atom_types[2] &&
         torsion_errors[ntorsion_errors].atom_types[2]==torsion_errors[i].atom_types[1] &&
         torsion_errors[ntorsion_errors].atom_types[3]==torsion_errors[i].atom_types[0] &&
         torsion_errors[ntorsion_errors].orig_env==torsion_errors[i].orig_env)) already_exists=TRUE;
      }

   /* IF IT DID NOT ALREADY EXIST, SAVE IT BY INCREMENTING THE NUMBER OF */
   /* KNOWN TORSION ERRORS (THIS PRESERVES THE NEW ERROR SINCE IT WAS ALREADY */
   /* ADDED TO END OF LIST PROVISIONALLY WHEN READ FROM THE .OUT FILE). */
   if (!already_exists) ntorsion_errors++;

   return;
   }

void save_if_new_bond_error()
/*============================================================================*/
/* PURPOSE: DETERMINE IT THE BOND ERROR JUST ADDED TO THE END OF THE LIST
/*	IS A NEW ERROR.  IF SO, SAVE THE ERROR BY INCREMENTING THE NUMBER OF
/*	KNOWN BOND ERRORS.
*/
   {
   int already_exists;	/* FLAG; T IF BOND ERROR ALREADY IN LIST. */
   int i;		/* LOOP INDEX; WHICH BOND ERROR. */

   already_exists=FALSE;
   for (i=0; i<nbond_errors; i++)
      {
      if ((bond_errors[nbond_errors].atom_types[0]==bond_errors[i].atom_types[0] &&
         bond_errors[nbond_errors].atom_types[1]==bond_errors[i].atom_types[1]) ||
         (bond_errors[nbond_errors].atom_types[0]==bond_errors[i].atom_types[1] &&
         bond_errors[nbond_errors].atom_types[1]==bond_errors[i].atom_types[0])) already_exists=TRUE;
      }

   /* IF IT DID NOT ALREADY EXIST, SAVE IT BY INCREMENTING THE NUMBER OF */
   /* KNOWN BOND ERRORS (THIS PRESERVES THE NEW ERROR SINCE IT WAS ALREADY */
   /* ADDED TO END OF LIST PROVISIONALLY WHEN READ FROM THE .OUT FILE). */
   if (!already_exists) nbond_errors++;

   return;
   }

void save_if_new_angle_error()
/*============================================================================*/
/* PURPOSE: DETERMINE IF THE ANGLE ERROR JUST ADDED TO THE END OF THE LIST
/*	IS A NEW ERROR.  IF SO, SAVE THE ERROR BY INCREMENTING THE NUMBER OF
/*	KNOWN ANGLE ERRORS.
*/
   {
   int already_exists;	/* FLAG; T IF ANGLE ERROR ALREADY IN LIST. */
   int i;		/* LOOP INDEX; WHICH ANGLE ERROR. */

   already_exists=FALSE;
   for (i=0; i<nangle_errors; i++)
      {
      if ((angle_errors[nangle_errors].atom_types[0]==angle_errors[i].atom_types[0] &&
         angle_errors[nangle_errors].atom_types[1]==angle_errors[i].atom_types[1] &&
         angle_errors[nangle_errors].atom_types[2]==angle_errors[i].atom_types[2] &&
         angle_errors[nangle_errors].orig_env==angle_errors[i].orig_env) ||
         (angle_errors[nangle_errors].atom_types[0]==angle_errors[i].atom_types[2] &&
         angle_errors[nangle_errors].atom_types[1]==angle_errors[i].atom_types[1] &&
         angle_errors[nangle_errors].atom_types[2]==angle_errors[i].atom_types[0] &&
         angle_errors[nangle_errors].orig_env==angle_errors[i].orig_env)) already_exists=TRUE;
      }

   /* IF IT DID NOT ALREADY EXIST, SAVE IT BY INCREMENTING THE NUMBER OF */
   /* KNOWN ANGLE ERRORS (THIS PRESERVES THE NEW ERROR SINCE IT WAS ALREADY */
   /* ADDED TO END OF LIST PROVISIONALLY WHEN READ FROM THE .OUT FILE). */
   if (!already_exists) nangle_errors++;

   return;
   }

void save_if_new_oop_error()
/*============================================================================*/
/* PURPOSE: DETERMINE IF THE OUT OF PLANE ERROR JUST ADDED TO THE END OF THE LIST
/*	IS A NEW ERROR.  IF SO, SAVE THE ERROR BY INCREMENTING THE NUMBER OF
/*	KNOWN OUT OF PLANE ERRORS.
*/
   {
   int already_exists;	/* FLAG; T IF OUT OF PLANE ERROR ALREADY IN LIST. */
   int i;		/* LOOP INDEX; WHICH OUT OF PLANE ERROR. */

   already_exists=FALSE;
   for (i=0; i<noutofplane_errors; i++)
      {
      if (outofplane_errors[noutofplane_errors].atom_types[0]==outofplane_errors[i].atom_types[0] &&
         outofplane_errors[noutofplane_errors].atom_types[1]==outofplane_errors[i].atom_types[1] &&
         outofplane_errors[noutofplane_errors].atom_types[2]==outofplane_errors[i].atom_types[2]) already_exists=TRUE;
      }

   /* IF IT DID NOT ALREADY EXIST, SAVE IT BY INCREMENTING THE NUMBER OF */
   /* KNOWN OOP ERRORS (THIS PRESERVES THE NEW ERROR SINCE IT WAS ALREADY */
   /* ADDED TO END OF LIST PROVISIONALLY WHEN READ FROM THE .OUT FILE). */
   if (!already_exists) noutofplane_errors++;

   return;
   }

int make_parameters()
/*============================================================================*/
/* PURPOSE: PROCESS EACH ERROR BY GENERATING THE MISSING PARAMETERS BASED ON
/*	SIMILARITIES TO EXISTING ATOM TYPE COMBINATIONS, OR ASK USER IF
/*	INSUFFICIENT SIMILARITIES EXIST.
*/
   {
   /* STOP PROCESSING IF USER ABORTS AT ANY POINT. */
   if (!process_torsion_errors()) return(FALSE);
   if (!process_bond_errors()) return(FALSE);
   if (!process_angle_errors()) return(FALSE);
   if (!process_outofplane_errors()) return(FALSE);

   /* SUCCESSFUL IF ANY ERRORS WERE PROCESSED. */
   if (ntorsion_errors || nbond_errors || nangle_errors || noutofplane_errors) return(TRUE);
   else return(FALSE);
   }

int process_torsion_errors()
/*============================================================================*/
/* PURPOSE: GENERATE MISSING TORSION PARAMETERS BASED ON SIMILARITIES WITH
/*	EXISTING ATOM TYPE COMBINATIONS, OR ASK USER IF INSUFFICIENT
/*	SIMILARITIES EXIST.
*/
   {
   int torsion[4];	/* ATOM TYPES IN TORSION. */
   int environment;	/* STRUCTURAL ENVIRONMENT OF TORSION. */
   double default_v1;	/* TORSIONAL V1 PARAMETER. */
   double default_v2;	/* TORSIONAL V2 PARAMETER. */
   double default_v3;	/* TORSIONAL V3 PARAMETER. */
   int match;		/* TYPE OF MATCH FOUND FOR THIS PARAMETER. */
   char prompt[80];	/* PROMPT STRING TO USER. */
   int i;		/* LOOP INDEX. */
   int j;		/* LOOP INDEX. */
   int status;		/* USER IO STATUS. */

   for (i=0; i<ntorsion_errors; i++)
      {
      /* DEFINE THE ORIGINAL TORSION, THEN LOOK FOR A SIMILAR TORSION. */
      for (j=0; j<4; j++)
         {
         torsion[j]=torsion_errors[i].atom_types[j];
         }
      environment=torsion_errors[i].orig_env;
      match=similar_torsion_exists(torsion,&environment,&default_v1,&default_v2,&default_v3);

      /* PARAMETER VALUE WILL BE SAVED ONLY IF *NOT* AN EXACT MATCH. */
      torsion_errors[i].save_param=(match==SIMILAR_MATCH || match==NO_MATCH);

      if (auto_take || match==EXACT_MATCH)
         {
         /* AUTOMATICALLY TAKE ALL DEFAULTS. */
         torsion_errors[i].v1=default_v1;
         torsion_errors[i].v2=default_v2;
         torsion_errors[i].v3=default_v3;
         }
      else
         {
         /* ALLOW USER TO EDIT ONLY NON IDENTICALLY MATCHED PARAMETERS. */
         /* TELL THE USER WHICH TORSION ANGLE WE ARE WORKING WITH. */
         printf("\nTorsion %d-%d-%d-%d",
            torsion_errors[i].atom_types[0],
            torsion_errors[i].atom_types[1],
            torsion_errors[i].atom_types[2],
            torsion_errors[i].atom_types[3]);
         if (torsion_errors[i].orig_env==FOUR_RING) printf(" (4 ring)");
         if (match==SIMILAR_MATCH)
            printf(":   (defaults from similar torsion %d-%d-%d-%d)\n",
               torsion[0],torsion[1],torsion[2],torsion[3]);
         else
            printf(":   (no similar torsion exists)\n");

         /* RUN EACH DEFAULT PARAMETER BY THE USER FOR APPROVAL OR CHANGE. */
         sprintf(prompt,"Enter one-fold torsional constant [%.3f]",default_v1);
         status=gdouble(prompt,NULL_OK,RANGE_CHECK,-20.0,20.0,&torsion_errors[i].v1);
         if (status==ABORT_ENTRY) return(FALSE);
         if (status==NULL_ENTRY) torsion_errors[i].v1=default_v1;

         sprintf(prompt,"Enter two-fold torsional constant [%.3f]",default_v2);
         status=gdouble(prompt,NULL_OK,RANGE_CHECK,-100.0,100.0,&torsion_errors[i].v2);
         if (status==ABORT_ENTRY) return(FALSE);
         if (status==NULL_ENTRY) torsion_errors[i].v2=default_v2;

         sprintf(prompt,"Enter three-fold torsional constant [%.3f]",default_v3);
         status=gdouble(prompt,NULL_OK,RANGE_CHECK,-20.0,20.0,&torsion_errors[i].v3);
         if (status==ABORT_ENTRY) return(FALSE);
         if (status==NULL_ENTRY) torsion_errors[i].v3=default_v3;
         }
      }

   return(TRUE);
   }

int get_torsion_environment(torsion_numbers)
/*============================================================================*/
/* PURPOSE: DETERMINE THE STRUCTURAL ENVIRONMENT OF THE SPECIFIED TORSION.
/*	BOTH THE CONNECT STRUCTURE AND THE SPECIFIED ATOM NUMBERS INDEX FROM 0.
*/
   int torsion_numbers[4];	/* ATOM NUMBERS OF TORSION. */
   {
   int i;			/* LOOP INDEX. */

   /* TORSION IS IN A 4-MEMBERED RING IFF ATOM1 IS BONDED TO ATOM4, OTHERWISE */
   /* IT IS IN A NORMAL ENVIRONMENT, SO CHECK EACH BOND FROM FIRST ATOM. */
   for (i=0; i<connect[torsion_numbers[0]].nbonds; i++)
      {
      if (connect[torsion_numbers[0]].bonds[i]==torsion_numbers[3]) return(FOUR_RING);
      }

   return(NORMAL);	/* WAS NOT IN A 4-MEMBERED RING, MUST BE NORMAL. */
   }

int similar_torsion_exists(torsion,environment,default_v1,default_v2,default_v3)
/*============================================================================*/
/* PURPOSE: DETERMINE THE TRANSFORMED TORSION MOST SIMILAR TO THE ORIGINAL
/*	TORSION.  IF SUCH A TRANSFORMED TORSION EXISTS, THE TORSIONAL CONSTANTS
/*	ARE RETRIEVED FROM THE SIMILAR DATA SET; IF NOT, THEY ARE SET TO DEFAULT
/*	VALUES.
*/
   int torsion[4];		/* ATOM TYPES IN TORSION. */
   int* environment;		/* STRUCTURAL ENVIRONMENT OF TORSION. */
   double* default_v1;		/* NEW TORSIONAL V1 PARAMETER. */
   double* default_v2;		/* NEW TORSIONAL V2 PARAMETER. */
   double* default_v3;		/* NEW TORSIONAL V3 PARAMETER. */
   {
   int orig_environment;	/* ORIGINAL ENVIRONMENT OF TORSION. */
   int new_torsion[4][4];	/* COPY OF TORSION WITH NTH ATOM TRANSFORMED. */
   int new_environment[4];	/* STRUCTURAL ENVIRONMENT OF NEW_TORSION[N]. */
   int valid_transform[4];	/* FLAG; TRANSFORM OF NTH ATOM IS DEFINED. */
   int cost[4];			/* COST OF ACHIEVING NEW_TORSION[N]. */
   int most_similar;		/* WHICH TRANSFORM IS MOST SIMILAR TO ORIG. */
   int i;			/* LOOP INDEX. */
   int j;			/* LOOP INDEX. */

   /* SET UP THE LAST DITCH DEFAULTS IN CASE NO SIMILAR TORSION EXISTS. */
   *default_v1=0.0;
   *default_v2=0.0;
   *default_v3=0.0;

   /* SINCE WE KNOW MORE PARAMETERS THAN MM2, CHECK FOR AN EXACT MATCH FIRST. */
   orig_environment= *environment;
   if (valid_torsion(torsion,environment,default_v1,default_v2,default_v3))
      {
      if (debug_on) printf("\n   Torsion %d-%d-%d-%d found exact match (env=%d/%d).\n",
         torsion[0],torsion[1],torsion[2],torsion[3],orig_environment,*environment);
      return(EXACT_MATCH);
      }

   /* NO EXACT MATCH, SO TRY TRANSFORMING ONE ATOM IN EACH OF FOUR COPIES OF */
   /* TORSION, THEN PICK THE NEW TORSION MOST SIMILAR TO THE ORIGINAL TORSION. */
   if (debug_on) printf("\n   Torsion %d-%d-%d-%d transformations:",torsion[0],torsion[1],torsion[2],torsion[3]);
   for (i=0; i<4; i++)
      {
      for (j=0; j<4; j++)
         {
         new_torsion[i][j]=torsion[j];
         }
      new_environment[i]= *environment;
      valid_transform[i]=find_valid_torsion(torsion,i,&new_environment[i],new_torsion[i],&cost[i]);
      }

   /* IF A NEW TORSION WHICH TRANSFORMED EITHER OF THE LEAST IMPORTANT ATOMS */
   /* (THESE ARE THE OUTSIDE ATOMS, NUMBERS 0 AND 3) EXISTS, WE CHOOSE THAT */
   /* TRANSFORMED TORSION MOST SIMILAR TO THE ORIGINAL TORSION (IE, WHICHEVER */
   /* HAS THE MINIMUM TRANSFORMATION COST).  NOTE THAT NEW_ENVIRONMENT MAY */
   /* CHANGE TO REFLECT WHERE PARAMETERS HAVE BEEN FOUND. */
   if (valid_transform[0] && valid_transform[3])
      {
      /* BOTH OUTSIDE ATOMS SUCCESSFULLY TRANSFORMED, PICK THE BEST. */
      if (cost[0] < cost[3]) most_similar=0;
      else if (cost[3] < cost[0]) most_similar=3;
      else
         {
         /* BOTH EXIST WITH SAME COST, PICK LOWEST ORIGINAL ATOM TYPE. */
         if (torsion[0] <= torsion[3]) most_similar=0;
         else most_similar=3;
         }
      }
   else if (valid_transform[0]) most_similar=0;	/* ONLY ATOM 0 TRANSFORMED. */
   else if (valid_transform[3]) most_similar=3;	/* ONLY ATOM 3 TRANSFORMED. */
   else
      {
      /* SINCE NEITHER OF THE OUTSIDE ATOMS WERE SUCCESSFULLY TRANSFORMED, */
      /* WE NOW LOOK AT TRANSFORMATIONS OF THE MORE IMPORTANT CENTER ATOMS */
      /* (THESE ARE THE INSIDE ATOMS, NUMBERS 1 AND 2). */
      if (valid_transform[1] && valid_transform[2])
         {
         if (cost[1] < cost[2]) most_similar=1;
         else if (cost[2] < cost[1]) most_similar=2;
         else
            {
            /* BOTH EXIST AND COST THE SAME, PICK LOWEST ORIGINAL ATOM TYPE. */
            if (torsion[1] <= torsion[2]) most_similar=1;
            else most_similar=2;
            }
         }
      else if (valid_transform[1]) most_similar=1;	/* ONLY ATOM 1 TRANSFORMED. */
      else if (valid_transform[2]) most_similar=2;	/* ONLY ATOM 2 TRANSFORMED. */
      else return(NO_MATCH);	/* NONE OF FOUR ATOMS TRANSFORMED! */
      }

   /* IF WE DROPPED THROUGH TO HERE, WE HAVE A VALID SIMILAR TORSION, SO */
   /* COPY THE SINGLE CHANGED ATOM TYPE INTO THE ORIGINAL TORSION AND */
   /* CALL VALID_TORSION SOLELY TO RETRIEVE THE DEFAULT PARAMETERS. */
   torsion[most_similar]=new_torsion[most_similar][most_similar];
   valid_torsion(torsion,&new_environment[most_similar],default_v1,default_v2,default_v3);
   *environment=new_environment[most_similar];
   return(SIMILAR_MATCH);
   }

int find_valid_torsion(orig_torsion,which_atom,environment,new_torsion,cost)
/*============================================================================*/
/* PURPOSE: REPEATEDLY TRANSFORM THE SELECTED ATOM OF THE ORIGINAL TORSION
/*	UNTIL WE EITHER FIND A VALID TORSION FOR WHICH PARAMETERS EXIST OR WE
/*	RUN OUT OF TRANSFORMATIONS.  NOTE THAT WE START LOOKING FOR A MATCH IN
/*	THE DATA FOR THE MATCHING STRUCTURAL ENVIRONMENT, BUT IF IT IS NOT
/*	FOUND THERE, WE ALWAYS DEFAULT TO LOOKING IN THE "NORMAL" ENVIRONMENT.
/*
*/
   int orig_torsion[4];		/* ORIGINAL TORSION ANGLE ATOM TYPES. */
   int which_atom;		/* WHICH ATOM TO TRANSFORM IN TORSION. */
   int* environment;		/* STRUCTURAL ENVIRONMENT OF TORSION. */
   int new_torsion[4];		/* NEW, TRANSFORMED TORSION IF SUCCESSFUL. */
   int* cost;			/* COST OF TRANSFORMATION. */
   {
   int orig_environment;	/* ORIGINAL DESIRED ENVIRONMENT. */
   int transform_level;		/* WHICH TRANSFORMATION LEVEL TO USE. */
   int valid_transform;		/* FLAG; WHETHER A VALID TRANSFORM WAS FOUND. */
   double v1;			/* V1 PARAMETER OF NEW TORSION (JUNK HERE). */
   double v2;			/* V2 PARAMETER OF NEW TORSION (JUNK HERE). */
   double v3;			/* V3 PARAMETER OF NEW TORSION (JUNK HERE). */

   /* KEEP TRANSFORMING THE SELECTED ATOM UNTIL WE EITHER FIND A VALID */
   /* TORSION FOR WHICH PARAMETERS EXIST OR WE RUN OUT OF TRANSFORMATIONS. */
   orig_environment= *environment;	/* REMEMBER ORIGINAL ENVIRONMENT. */
   transform_level=0;			/* ALWAYS START AT THE BEGINNING. */
   while ((valid_transform=transform(orig_torsion[which_atom],&transform_level,&new_torsion[which_atom],cost)) &&
      !valid_torsion(new_torsion,environment,&v1,&v2,&v3))
      {
      if (debug_on)
         {
         printf("\n   Transforming atom #%d yields %2d-%2d-%2d-%2d: no parameters (env=%d/?).",
            which_atom+1,new_torsion[0],new_torsion[1],new_torsion[2],new_torsion[3],orig_environment);
         }
      }
   if (*environment!=orig_environment) (*cost)+=5;	/* ADDITIONAL COST. */

   if (debug_on)
      {
      if (valid_transform) printf("\n   Transforming atom #%d yields %2d-%2d-%2d-%2d: parameters exist (cost=%d,env=%d/%d).",
         which_atom+1,new_torsion[0],new_torsion[1],new_torsion[2],new_torsion[3],*cost,orig_environment,*environment);
      printf("\n");
      }

   return(valid_transform);
   }

int valid_torsion(torsion,environment,v1,v2,v3)
/*============================================================================*/
/* PURPOSE: LOOK UP THE SPECIFIED TORSION ANGLE IN THE LIST OF KNOWN TORSIONS
/*	FOR THE ENVIRONMENT; IF IT EXISTS, RETURN THE ASSOCIATED PARAMETERS.
*/
   int torsion[4];	/* TORSION ANGLE TO LOOK UP. */
   int* environment;	/* TORSION ENVIRONMENT (4 MEMBERED RING, ETC). */
   double* v1;		/* V1 PARAMETER OF MATCHED TORSION, IF ANY. */
   double* v2;		/* V2 PARAMETER OF MATCHED TORSION, IF ANY. */
   double* v3;		/* V3 PARAMETER OF MATCHED TORSION, IF ANY. */
   {
   int i;		/* LOOP INDEX. */

   if (*environment==FOUR_RING)
      {
      /* LOOK FOR THE TORSION IN THE KNOWN 4-MEMBERED RING TORSIONS. */
      for (i=0; i<ntorsion_params_4ring; i++)
         {
         if ((torsion_params_4ring[i].atom_types[0]==torsion[0] &&
            torsion_params_4ring[i].atom_types[1]==torsion[1] &&
            torsion_params_4ring[i].atom_types[2]==torsion[2] &&
            torsion_params_4ring[i].atom_types[3]==torsion[3]) ||
            (torsion_params_4ring[i].atom_types[0]==torsion[3] &&
            torsion_params_4ring[i].atom_types[1]==torsion[2] &&
            torsion_params_4ring[i].atom_types[2]==torsion[1] &&
            torsion_params_4ring[i].atom_types[3]==torsion[0]))
            {
            /* FOUND THE SPECIFIED TORSION, COPY PARAMETERS. */
            *v1=torsion_params_4ring[i].v1;
            *v2=torsion_params_4ring[i].v2;
            *v3=torsion_params_4ring[i].v3;
            return(TRUE);
            }
         }
      }

   if (*environment==NORMAL || *environment==FOUR_RING)
      {
      /* LOOK FOR THE TORSION IN THE KNOWN NORMAL TORSIONS. */
      for (i=0; i<ntorsion_params_normal; i++)
         {
         if ((torsion_params_normal[i].atom_types[0]==torsion[0] &&
            torsion_params_normal[i].atom_types[1]==torsion[1] &&
            torsion_params_normal[i].atom_types[2]==torsion[2] &&
            torsion_params_normal[i].atom_types[3]==torsion[3]) ||
            (torsion_params_normal[i].atom_types[0]==torsion[3] &&
            torsion_params_normal[i].atom_types[1]==torsion[2] &&
            torsion_params_normal[i].atom_types[2]==torsion[1] &&
            torsion_params_normal[i].atom_types[3]==torsion[0]))
            {
            /* FOUND THE SPECIFIED TORSION, COPY PARAMETERS. */
            *v1=torsion_params_normal[i].v1;
            *v2=torsion_params_normal[i].v2;
            *v3=torsion_params_normal[i].v3;
            *environment=NORMAL;	/* MAY HAVE CHANGED ENVIRONMENT. */
            return(TRUE);
            }
         }
      }

   /* NO MATCH FOUND IN THE SPECIFIED ENVIRONMENT. */
   return(FALSE);
   }

int process_bond_errors()
/*============================================================================*/
/* PURPOSE: GENERATE MISSING BOND PARAMETERS BASED ON SIMILARITIES WITH
/*	EXISTING ATOM TYPE COMBINATIONS, OR ASK USER IF INSUFFICIENT
/*	SIMILARITIES EXIST.
*/
   {
   int bond[2];		/* ATOM TYPES IN BOND. */
   double default_s;	/* BOND STRETCHING PARAMETER. */
   double default_l0;	/* BOND MINIMUM ENERGY LENGTH PARAMETER. */
   double default_bmom;	/* BOND MOMENT PARAMETER. */
   double default_slps;	/* BOND SLOPE STRETCHING EQUATION PARAMETER. */
   double default_slpt;	/* BOND SLOPE LENGTH EQUATION PARAMETER. */
   int match;		/* TYPE OF MATCH FOUND FOR THIS PARAMETER. */
   char prompt[80];	/* PROMPT STRING TO USER. */
   int i;		/* LOOP INDEX. */
   int status;		/* USER IO STATUS. */

   for (i=0; i<nbond_errors; i++)
      {
      /* DEFINE THE ORIGINAL BOND, THEN LOOK FOR A SIMILAR BOND. */
      bond[0]=bond_errors[i].atom_types[0];
      bond[1]=bond_errors[i].atom_types[1];
      match=similar_bond_exists(bond,&default_s,&default_l0,&default_bmom,&default_slps,&default_slpt);

      /* PARAMETER VALUE WILL BE SAVED ONLY IF *NOT* AN EXACT MATCH. */
      bond_errors[i].save_param=(match==SIMILAR_MATCH || match==NO_MATCH);

      if (auto_take || match==EXACT_MATCH)
         {
         /* AUTOMATICALLY TAKE ALL DEFAULTS. */
         bond_errors[i].s=default_s;
         bond_errors[i].l0=default_l0;
         bond_errors[i].bmom=default_bmom;
         bond_errors[i].slps=default_slps;
         bond_errors[i].slpt=default_slpt;
         }
      else
         {
         /* ALLOW USER TO EDIT ONLY NON IDENTICALLY MATCHED PARAMETERS. */
         /* TELL THE USER WHICH BOND WE ARE WORKING WITH. */
         printf("\nBond %d-%d",
            bond_errors[i].atom_types[0],
            bond_errors[i].atom_types[1]);
         if (match==SIMILAR_MATCH)
            printf(":   (defaults from similar bond %d-%d)\n",bond[0],bond[1]);
         else
            printf(":   (no similar bond exists)\n");

         /* RUN EACH DEFAULT PARAMETER BY THE USER FOR APPROVAL OR CHANGE. */
         sprintf(prompt,"Enter stretching constant [%.3f]",default_s);
         status=gdouble(prompt,NULL_OK,RANGE_CHECK,0.0,30.0,&bond_errors[i].s);
         if (status==ABORT_ENTRY) return(FALSE);
         if (status==NULL_ENTRY) bond_errors[i].s=default_s;

         sprintf(prompt,"Enter minimum energy bond length [%.3f]",default_l0);
         status=gdouble(prompt,NULL_OK,RANGE_CHECK,0.5,3.0,&bond_errors[i].l0);
         if (status==ABORT_ENTRY) return(FALSE);
         if (status==NULL_ENTRY) bond_errors[i].l0=default_l0;

         sprintf(prompt,"Enter bond moment constant [%.3f]",default_bmom);
         status=gdouble(prompt,NULL_OK,RANGE_CHECK,-5.0,5.0,&bond_errors[i].bmom);
         if (status==ABORT_ENTRY) return(FALSE);
         if (status==NULL_ENTRY) bond_errors[i].bmom=default_bmom;

         sprintf(prompt,"Enter slope of bond order stretching equation [%.3f]",default_slps);
         status=gdouble(prompt,NULL_OK,NO_RANGE_CHECK,0.0,0.0,&bond_errors[i].slps);
         if (status==ABORT_ENTRY) return(FALSE);
         if (status==NULL_ENTRY) bond_errors[i].slps=default_slps;

         sprintf(prompt,"Enter slope of bond order length equation [%.3f]",default_slpt);
         status=gdouble(prompt,NULL_OK,NO_RANGE_CHECK,0.0,0.0,&bond_errors[i].slpt);
         if (status==ABORT_ENTRY) return(FALSE);
         if (status==NULL_ENTRY) bond_errors[i].slpt=default_slpt;
         }
      }

   return(TRUE);
   }

int similar_bond_exists(bond,default_s,default_l0,default_bmom,default_slps,default_slpt)
/*============================================================================*/
/* PURPOSE: DETERMINE IF A SIMILAR BOND EXISTS.  IF SO, THE
/*	BOND CONSTANTS ARE RETRIEVED FROM THE SIMILAR DATA SET; IF NOT,
/*	THE BOND CONSTANTS ARE SET TO DEFAULT VALUES.
*/
   int bond[2];			/* ATOM TYPES IN BOND. */
   double* default_s;		/* BOND STRETCHING PARAMETER. */
   double* default_l0;		/* BOND MINIMUM ENERGY LENGTH PARAMETER. */
   double* default_bmom;	/* BOND MOMENT PARAMETER. */
   double* default_slps;	/* BOND SLOPE STRETCHING EQUATION PARAMETER. */
   double* default_slpt;	/* BOND SLOPE LENGTH EQUATION PARAMETER. */
   {
   int new_bond[2][2];		/* COPY OF BOND WITH NTH ATOM TRANSFORMED. */
   int valid_transform[2];	/* FLAG; TRANSFORM OF NTH ATOM IS DEFINED. */
   int cost[2];			/* COST OF ACHIEVING NEW_BOND[N]. */
   int most_similar;		/* WHICH TRANSFORM IS MOST SIMILAR TO ORIG. */
   int i;			/* LOOP INDEX. */
   int j;			/* LOOP INDEX. */

   void define_bond_defaults();	/* DETERMINE LAST DITCH DEFAULTS. */

   /* SET UP THE LAST DITCH DEFAULTS IN CASE NO SIMILAR BOND EXISTS. */
   define_bond_defaults(bond,default_s,default_l0,default_bmom,default_slps,default_slpt);

   /* SINCE WE KNOW MORE PARAMETERS THAN MM2, CHECK FOR AN EXACT MATCH FIRST. */
   if (valid_bond(bond,default_s,default_l0,default_bmom,default_slps,default_slpt))
      {
      if (debug_on) printf("\n   Bond %d-%d found exact match.\n",bond[0],bond[1]);
      return(EXACT_MATCH);
      }

   /* NO EXACT MATCH, SO TRY TRANSFORMING ONE ATOM IN EACH OF TWO COPIES OF */
   /* BOND, THEN PICK THE NEW BOND MOST SIMILAR TO THE ORIGINAL BOND. */
   if (debug_on) printf("\n   Bond %d-%d transformations:",bond[0],bond[1]);
   for (i=0; i<2; i++)
      {
      for (j=0; j<2; j++)
         {
         new_bond[i][j]=bond[j];
         }
      valid_transform[i]=find_valid_bond(bond,i,new_bond[i],&cost[i]);
      }

   /* IF BOTH ATOMS SUCCESSFULLY TRANSFORMED, PICK THE BEST. */
   if (valid_transform[0] && valid_transform[1])
      {
      if (cost[0] < cost[1]) most_similar=0;
      else if (cost[1] < cost[0]) most_similar=1;
      else
         {
         /* BOTH EXIST WITH SAME COST, PICK LOWEST ORIGINAL ATOM TYPE. */
         if (bond[0] <= bond[1]) most_similar=0;
         else most_similar=1;
         }
      }
   else if (valid_transform[0]) most_similar=0;	/* ONLY ATOM 0 TRANSFORMED. */
   else if (valid_transform[1]) most_similar=1;	/* ONLY ATOM 1 TRANSFORMED. */
   else return(NO_MATCH);	/* NEITHER ATOM TRANSFORMED! */

   /* IF WE DROPPED THROUGH TO HERE, WE HAVE A VALID SIMILAR BOND, SO */
   /* COPY THE SINGLE CHANGED ATOM TYPE INTO THE ORIGINAL BOND AND */
   /* CALL VALID_BOND SOLELY TO RETRIEVE THE DEFAULT PARAMETERS. */
   bond[most_similar]=new_bond[most_similar][most_similar];
   valid_bond(bond,default_s,default_l0,default_bmom,default_slps,default_slpt);
   return(SIMILAR_MATCH);
   }

void define_bond_defaults(bond,default_s,default_l0,default_bmom,default_slps,default_slpt)
/*============================================================================*/
/* PURPOSE: DEFINE THE LAST DITCH DEFAULTS TO BE USED AS PARAMETERS FOR THIS
/*	BOND IF NO SIMILAR BOND EXISTS.
*/
   int bond[2];			/* ATOM TYPES IN BOND. */
   double* default_s;		/* BOND STRETCHING PARAMETER. */
   double* default_l0;		/* BOND MINIMUM ENERGY LENGTH PARAMETER. */
   double* default_bmom;	/* BOND MOMENT PARAMETER. */
   double* default_slps;	/* BOND SLOPE STRETCHING EQUATION PARAMETER. */
   double* default_slpt;	/* BOND SLOPE LENGTH EQUATION PARAMETER. */
   {
   *default_bmom=0.0;
   *default_slps=0.0;
   *default_slpt=0.0;

   if (bond_contains(bond,20))
      {
      *default_s=4.6;
      *default_l0=0.6;
      }
   else if (bond_contains(bond,5) || bond_contains(bond,21) || bond_contains(bond,23) ||
      bond_contains(bond,24) || bond_contains(bond,28) || bond_contains(bond,48))
      {
      *default_s=4.6;
      *default_l0=1.0;
      }
   else if (bond_contains(bond,2) || bond_contains(bond,3) || bond_contains(bond,9) ||
      bond_contains(bond,17) || bond_contains(bond,22) || bond_contains(bond,26) ||
      bond_contains(bond,29) || bond_contains(bond,37) || bond_contains(bond,38) ||
      bond_contains(bond,40) || bond_contains(bond,41) || bond_contains(bond,42) ||
      bond_contains(bond,46) || bond_contains(bond,49) || bond_contains(bond,50))
      {
      *default_s=9.6;
      *default_l0=1.35;
      }
   else if (bond_contains(bond,4) || bond_contains(bond,10) || bond_contains(bond,43) ||
      bond_contains(bond,45))
      {
      *default_s=15.0;
      *default_l0=1.2;
      }
   else
      {
      *default_s=4.4;
      *default_l0=1.49;
      }

   return;
   }

int bond_contains(bond,atom_type)
/*============================================================================*/
/* PURPOSE: DETERMINE IF THE BOND CONTAINS THE SPECIFIED ATOM TYPE IN EITHER
/*	POSITION.
*/
   int bond[2];		/* ATOM TYPES IN BOND. */
   int atom_type;	/* ATOM TYPE TO LOOK FOR. */
   {
   if (bond[0]==atom_type || bond[1]==atom_type) return(TRUE);
   else return(FALSE);
   }

int find_valid_bond(orig_bond,which_atom,new_bond,cost)
/*============================================================================*/
/* PURPOSE: REPEATEDLY TRANSFORM THE SELECTED ATOM OF THE ORIGINAL BOND
/*	UNTIL WE EITHER FIND A VALID BOND FOR WHICH PARAMETERS EXIST OR WE
/*	RUN OUT OF TRANSFORMATIONS.
*/
   int orig_bond[2];	/* ORIGINAL BOND ATOM TYPES. */
   int which_atom;	/* WHICH ATOM TO TRANSFORM IN BOND. */
   int new_bond[2];	/* NEW, TRANSFORMED BOND IF SUCCESSFUL. */
   int* cost;		/* COST OF TRANSFORMATION. */
   {
   int transform_level;	/* WHICH TRANSFORMATION LEVEL TO USE. */
   int valid_transform;	/* FLAG; WHETHER A VALID TRANSFORM WAS FOUND. */
   double s;		/* S PARAMETER OF NEW BOND (JUNK HERE). */
   double l0;		/* L0 PARAMETER OF NEW BOND (JUNK HERE). */
   double bmom;		/* BMOM PARAMETER OF NEW BOND (JUNK HERE). */
   double slps;		/* SLPS PARAMETER OF NEW BOND (JUNK HERE). */
   double slpt;		/* SLPT PARAMETER OF NEW BOND (JUNK HERE). */

   /* KEEP TRANSFORMING THE SELECTED ATOM UNTIL WE EITHER FIND A VALID */
   /* BOND FOR WHICH PARAMETERS EXIST OR WE RUN OUT OF TRANSFORMATIONS. */
   transform_level=0;	/* ALWAYS START AT THE BEGINNING. */
   while ((valid_transform=transform(orig_bond[which_atom],&transform_level,&new_bond[which_atom],cost)) &&
      !valid_bond(new_bond,&s,&l0,&bmom,&slps,&slpt))
      {
      if (debug_on)
         {
         printf("\n   Transforming atom #%d yields %2d-%2d: no parameters.",
            which_atom+1,new_bond[0],new_bond[1]);
         }
      }

   if (debug_on)
      {
      if (valid_transform) printf("\n   Transforming atom #%d yields %2d-%2d: parameters exist (cost=%d).",
         which_atom+1,new_bond[0],new_bond[1],*cost);
      printf("\n");
      }

   return(valid_transform);
   }

int valid_bond(bond,s,l0,bmom,slps,slpt)
/*============================================================================*/
/* PURPOSE: LOOK UP THE SPECIFIED BOND IN THE LIST OF KNOWN
/*	BONDS; IF IT EXISTS, RETURN THE ASSOCIATED PARAMETERS.
*/
   int bond[2];		/* BOND TO LOOK UP. */
   double* s;		/* S PARAMETER OF MATCHED BOND, IF ANY. */
   double* l0;		/* L0 PARAMETER OF MATCHED BOND, IF ANY. */
   double* bmom;	/* BMOM PARAMETER OF MATCHED BOND, IF ANY. */
   double* slps;	/* SLPS PARAMETER OF MATCHED BOND, IF ANY. */
   double* slpt;	/* SLPT PARAMETER OF MATCHED BOND, IF ANY. */
   {
   int i;		/* LOOP INDEX. */

   /* LOOK FOR THE SPECIFIED BOND IN THE LIST OF KNOWN BONDS. */
   for (i=0; i<nbond_params; i++)
      {
      if ((bond_params[i].atom_types[0]==bond[0] &&
         bond_params[i].atom_types[1]==bond[1]) ||
         (bond_params[i].atom_types[0]==bond[1] &&
         bond_params[i].atom_types[1]==bond[0]))
         {
         /* FOUND THE SPECIFIED BOND, COPY PARAMETERS. */
         *s=bond_params[i].s;
         *l0=bond_params[i].l0;
         *bmom=bond_params[i].bmom;
         *slps=bond_params[i].slps;
         *slpt=bond_params[i].slpt;
         return(TRUE);
         }
      }

   /* LOOKED THROUGH ALL KNOWN BONDS, NO MATCH FOUND. */
   return(FALSE);
   }

int process_angle_errors()
/*============================================================================*/
/* PURPOSE: GENERATE MISSING ANGLE PARAMETERS BASED ON SIMILARITIES WITH
/*	EXISTING ATOM TYPE COMBINATIONS, OR ASK USER IF INSUFFICIENT
/*	SIMILARITIES EXIST.
*/
   {
   int angle[3];	/* ATOM TYPES IN 2D ANGLE. */
   int environment;	/* STRUCTURAL ENVIRONMENT OF ANGLE. */
   double default_b;	/* BENDING CONSTANT PARAMETER. */
   double default_t;	/* MINIMUM ENERGY BOND ANGLE PARAMETER. */
   int match;		/* TYPE OF MATCH FOUND FOR THIS PARAMETER. */
   char prompt[80];	/* PROMPT STRING TO USER. */
   int i;		/* LOOP INDEX. */
   int j;		/* LOOP INDEX. */
   int status;		/* USER IO STATUS. */

   for (i=0; i<nangle_errors; i++)
      {
      /* DEFINE THE ORIGINAL ANGLE, AND LOOK FOR A SIMILAR ANGLE. */
      for (j=0; j<3; j++)
         {
         angle[j]=angle_errors[i].atom_types[j];
         }
      environment=angle_errors[i].orig_env;
      match=similar_angle_exists(angle,&environment,&default_b,&default_t);

      /* PARAMETER VALUE WILL BE SAVED ONLY IF *NOT* AN EXACT MATCH. */
      angle_errors[i].save_param=(match==SIMILAR_MATCH || match==NO_MATCH);

      if (auto_take || match==EXACT_MATCH)
         {
         /* AUTOMATICALLY TAKE ALL DEFAULTS. */
         angle_errors[i].b=default_b;
         angle_errors[i].t=default_t;
         }
      else
         {
         /* ALLOW USER TO EDIT ONLY NON IDENTICALLY MATCHED PARAMETERS. */
         /* TELL THE USER WHICH ANGLE WE ARE WORKING WITH. */
         printf("\nAngle %d-%d-%d",
            angle_errors[i].atom_types[0],
            angle_errors[i].atom_types[1],
            angle_errors[i].atom_types[2]);
         if (angle_errors[i].orig_env==THREE_RING) printf(" (3 ring)");
         else if (angle_errors[i].orig_env==FOUR_RING) printf(" (4 ring)");
         if (match==SIMILAR_MATCH)
            printf(":   (defaults from similar angle %d-%d-%d)\n",angle[0],angle[1],angle[2]);
         else
            printf(":   (no similar angle exists)\n");

         /* RUN EACH DEFAULT PARAMETER BY THE USER FOR APPROVAL OR CHANGE. */
         sprintf(prompt,"Enter bending constant [%.3f]",default_b);
         status=gdouble(prompt,NULL_OK,RANGE_CHECK,0.0,10.0,&angle_errors[i].b);
         if (status==ABORT_ENTRY) return(FALSE);
         if (status==NULL_ENTRY) angle_errors[i].b=default_b;

         sprintf(prompt,"Enter minimum energy bond angle [%.3f]",default_t);
         status=gdouble(prompt,NULL_OK,RANGE_CHECK,0.0,180.0,&angle_errors[i].t);
         if (status==ABORT_ENTRY) return(FALSE);
         if (status==NULL_ENTRY) angle_errors[i].t=default_t;
         }
      }

   return(TRUE);
   }

int get_angle_environment(angle_numbers,angle_types)
/*============================================================================*/
/* PURPOSE: DETERMINE THE STRUCTURAL ENVIRONMENT OF THE SPECIFIED ANGLE.
/*	ALL ATOM NUMBERS INDEX FROM ZERO, NOT ONE.
*/
   int angle_numbers[3];	/* ATOM NUMBERS OF ANGLE. */
   int angle_types[3];		/* ATOM TYPES OF ANGLE. */
   {
   int in_3_ring;		/* FLAG; IF IN 3-MEMBERED RING. */
   int in_4_ring;		/* FLAG; IF IN 4-MEMBERED RING. */
   int environment;		/* STRUCTURAL ENVIRONMENT. */
   double b;			/* B PARAMETER; JUNK HERE. */
   double t;			/* T PARAMETER; JUNK HERE. */
   int i;			/* LOOP INDEX. */
   int j;			/* LOOP INDEX. */

   /* WE HAVE A SLIGHTLY COMPLICATED SITUATION HERE.  ANGLE CAN BE EITHER */
   /* IN A 3 RING, A 4 RING, BOTH A 3 AND 4 RING, OR NEITHER.  SINCE THE */
   /* SEARCH FOR PARAMETERS IS BIFURCATED (3 RING, THEN NORMAL; *OR* 4 RING */
   /* THEN NORMAL) WITHOUT CROSSOVER BETWEEN 3 RING AND 4 RING, WE HAVE TO */
   /* DETERMINE RING MEMBERSHIPS FIRST, THEN, IF A MEMBER OF *BOTH* 3 AND 4 */
   /* RING CHECK IF THE 3 RING PARAMETERS EXIST (THEREFORE CAUSING US TO */
   /* RETURN "THREE_RING" AS THE ENVIRONMENT), OTHERWISE WE USE "FOUR_RING". */

   /* EXAMINE ALL THE BONDS FROM THE FIRST ATOM. */
   in_3_ring=in_4_ring=FALSE;
   for (i=0; i<connect[angle_numbers[0]].nbonds; i++)
      {
      /* IF FIRST BONDED TO THIRD, MUST BE IN A 3-MEMBERED RING. */
      if (connect[angle_numbers[0]].bonds[i]==angle_numbers[2]) in_3_ring=TRUE;

      /* CHECK BONDS OF THIRD; ANY CONNECTION WITH A (NON ANGLE MEMBER) ATOM */
      /* COMMON TO FIRST IMPLIES WE ARE IN A 4-MEMBERED RING. */
      for (j=0; j<connect[angle_numbers[2]].nbonds; j++)
         {
         if (connect[angle_numbers[2]].bonds[j]==connect[angle_numbers[0]].bonds[i] &&
            connect[angle_numbers[2]].bonds[j]!=angle_numbers[1]) in_4_ring=TRUE;
         }
      }

   /* NOW DECIDE WHAT TO DO WITH RAW ENVIRONMENT INFORMATION. */
   if (in_3_ring && in_4_ring)
      {
      environment=THREE_RING;
      if (!valid_angle(angle_types,&environment,&b,&t)) environment=FOUR_RING;
      }
   else if (in_3_ring) environment=THREE_RING;
   else if (in_4_ring) environment=FOUR_RING;
   else environment=NORMAL;

   return(environment);
   }

int similar_angle_exists(angle,environment,default_b,default_t)
/*============================================================================*/
/* PURPOSE: DETERMINE IF A TORSION ANGLE EXISTS.  IF SO, THE
/*	ANGLE CONSTANTS ARE RETRIEVED FROM THE SIMILAR DATA SET; IF NOT,
/*	THE ANGLE CONSTANTS ARE SET TO DEFAULT VALUES.
*/
   int angle[3];		/* ATOM TYPES IN 2D ANGLE. */
   int* environment;		/* STRUCTURAL ENVIRONMENT OF ANGLE. */
   double* default_b;		/* BENDING CONSTANT PARAMETER. */
   double* default_t;		/* MINIMUM ENERGY BOND ANGLE PARAMETER. */
   {
   int orig_environment;	/* ORIGINAL ENVIRONMENT OF TORSION. */
   int new_angle[3][3];		/* ANGLE WITH NTH ATOM TRANSFORMED. */
   int new_environment[3];	/* STRUCTURAL ENVIRONMENT OF NEW_ANGLE[N]. */
   int valid_transform[3];	/* FLAG; TRANSFORM OF NTH ATOM IS DEFINED. */
   int cost[3];			/* COST OF ACHIEVING NEW_ANGLE[N]. */
   int most_similar;		/* WHICH TRANSFORM IS MOST SIMILAR TO ORIG. */
   int i;			/* LOOP INDEX. */
   int j;			/* LOOP INDEX. */

   void define_angle_defaults();	/* DEFINE LAST DITCH DEFAULTS. */

   /* SET UP THE LAST DITCH DEFAULTS IN CASE NO SIMILAR ANGLE EXISTS. */
   define_angle_defaults(angle,default_b,default_t);

   /* SINCE WE KNOW MORE PARAMETERS THAN MM2, CHECK FOR AN EXACT MATCH FIRST. */
   orig_environment= *environment;
   if (valid_angle(angle,environment,default_b,default_t))
      {
      if (debug_on) printf("\n   Angle %d-%d-%d found exact match (env=%d/%d).\n",
         angle[0],angle[1],angle[2],orig_environment,*environment);
      return(EXACT_MATCH);
      }

   /* NO EXACT MATCH, SO TRY TRANSFORMING ONE ATOM IN EACH OF THREE COPIES OF */
   /* ANGLE, THEN PICK THE NEW ANGLE MOST SIMILAR TO THE ORIGINAL ANGLE. */
   if (debug_on) printf("\n   Angle %d-%d-%d transformations:",angle[0],angle[1],angle[2]);
   for (i=0; i<3; i++)
      {
      for (j=0; j<3; j++)
         {
         new_angle[i][j]=angle[j];
         }
      new_environment[i]= *environment;
      valid_transform[i]=find_valid_angle(angle,i,&new_environment[i],new_angle[i],&cost[i]);
      }

   /* IF A NEW ANGLE WHICH TRANSFORMED EITHER OF THE LEAST IMPORTANT ATOMS */
   /* (THESE ARE THE OUTSIDE ATOMS, NUMBERS 0 AND 2) EXISTS, WE CHOOSE THAT */
   /* TRANSFORMED ANGLE MOST SIMILAR TO THE ORIGINAL ANGLE (IE, WHICHEVER */
   /* HAS THE MINIMUM TRANSFORMATION COST).  NOTE THAT NEW_ENVIRONMENT MAY */
   /* CHANGE TO REFLECT WHERE PARAMETERS HAVE BEEN FOUND. */
   if (valid_transform[0] && valid_transform[2])
      {
      if (cost[0] < cost[2]) most_similar=0;
      else if (cost[2] < cost[0]) most_similar=2;
      else
         {
         /* BOTH EXIST AND COST THE SAME, PICK LOWEST ORIGINAL ATOM TYPE. */
         if (angle[0] <= angle[2]) most_similar=0;
         else most_similar=2;
         }
      }
   else if (valid_transform[0]) most_similar=0;	/* ONLY ATOM 0 TRANSFORMED. */
   else if (valid_transform[2]) most_similar=2;	/* ONLY ATOM 2 TRANSFORMED. */
   else
      {
      /* NEITHER OF THE OUTSIDE ATOMS WERE SUCCESSFULLY TRANSFORMED, */
      /* SO NOW WE LOOK AT TRANSFORMATIONS OF THE SINGLE CENTER ATOM. */
      if (valid_transform[1]) most_similar=1;
      else return(NO_MATCH);	/* NONE OF THREE ATOMS TRANSFORMED! */
      }

   /* IF WE DROPPED THROUGH TO HERE, WE HAVE A VALID SIMILAR ANGLE, SO */
   /* COPY THE SINGLE CHANGED ATOM TYPE INTO THE ORIGINAL ANGLE AND */
   /* CALL VALID_ANGLE SOLELY TO RETRIEVE THE DEFAULT PARAMETERS. */
   angle[most_similar]=new_angle[most_similar][most_similar];
   valid_angle(angle,&new_environment[most_similar],default_b,default_t);
   *environment=new_environment[most_similar];
   return(SIMILAR_MATCH);
   }

void define_angle_defaults(angle,default_b,default_t)
/*============================================================================*/
/* PURPOSE: DEFINE THE LAST DITCH DEFAULTS TO BE USED AS PARAMETERS FOR THIS
/*	ANGLE IF NO SIMILAR ANGLE EXISTS.
*/
   int angle[3];		/* ATOM TYPES IN 2D ANGLE. */
   double* default_b;		/* BENDING CONSTANT PARAMETER. */
   double* default_t;		/* MINIMUM ENERGY BOND ANGLE PARAMETER. */
   {
   if (angle[1]==1 || angle[1]==6 || angle[1]==8 || angle[1]==18 || angle[1]==19 || angle[1]==39)
      {
      *default_b=0.4;
      *default_t=107.0;
      }
   else if (angle[1]==15 || angle[1]==16 ||angle[1]==25)
      {
      *default_b=0.4;
      *default_t=96.0;
      }
   else if (angle[1]==2 || angle[1]==3 || angle[1]==9 || angle[1]==17 ||
      angle[1]==22 || angle[1]==26 || angle[1]==29 || angle[1]==37 ||
      angle[1]==38 || angle[1]==40 || angle[1]==41 || angle[1]==42 ||
      angle[1]==46 || angle[1]==49 || angle[1]==50)
      {
      *default_b=0.45;
      *default_t=120.0;
      }
   else if (angle[1]==4 || angle[1]==10 || angle[1]==43 || angle[1]==45)
      {
      *default_b=0.325;
      *default_t=180.0;
      }
   else
      {
      /* MUST HAVE BEEN AN INVALID INPUT ATOM TYPE - THIS IS A FATAL ERROR. */
      fprintf(stderr,"Error: Invalid atom type (%d) found as central angle atom.\n",angle[1]);
      exit(PROGRAM_FAILED);
      }

   return;
   }

int find_valid_angle(orig_angle,which_atom,environment,new_angle,cost)
/*============================================================================*/
/* PURPOSE: REPEATEDLY TRANSFORM THE SELECTED ATOM OF THE ORIGINAL ANGLE
/*	UNTIL WE EITHER FIND A VALID ANGLE FOR WHICH PARAMETERS EXIST OR WE
/*	RUN OUT OF TRANSFORMATIONS.
*/
   int orig_angle[3];		/* ORIGINAL ANGLE ATOM TYPES. */
   int which_atom;		/* WHICH ATOM TO TRANSFORM IN ANGLE. */
   int* environment;		/* STRUCTURAL ENVIRONMENT OF ANGLE. */
   int new_angle[3];		/* NEW, TRANSFORMED ANGLE IF SUCCESSFUL. */
   int* cost;			/* COST OF TRANSFORMATION. */
   {
   int orig_environment;	/* ORIGINAL DESIRED ENVIRONMENT. */
   int transform_level;		/* WHICH TRANSFORMATION LEVEL TO USE. */
   int valid_transform;		/* FLAG; WHETHER A VALID TRANSFORM WAS FOUND. */
   double b;			/* B PARAMETER OF NEW ANGLE (JUNK HERE). */
   double t;			/* T PARAMETER OF NEW ANGLE (JUNK HERE). */

   /* KEEP TRANSFORMING THE SELECTED ATOM UNTIL WE EITHER FIND A VALID */
   /* ANGLE FOR WHICH PARAMETERS EXIST OR WE RUN OUT OF TRANSFORMATIONS. */
   orig_environment= *environment;	/* REMEMBER ORIGINAL ENVIRONMENT. */
   transform_level=0;			/* ALWAYS START AT THE BEGINNING. */
   while ((valid_transform=transform(orig_angle[which_atom],&transform_level,&new_angle[which_atom],cost)) &&
      !valid_angle(new_angle,environment,&b,&t))
      {
      if (debug_on)
         {
         printf("\n   Transforming atom #%d yields %2d-%2d-%2d: no parameters (env=%d/?).",
            which_atom+1,new_angle[0],new_angle[1],new_angle[2],orig_environment);
         }
      }
   if (*environment!=orig_environment) (*cost)+=5;	/* ADDITIONAL COST. */

   if (debug_on)
      {
      if (valid_transform) printf("\n   Transforming atom #%d yields %2d-%2d-%2d: parameters exist (cost=%d,env=%d/%d).",
         which_atom+1,new_angle[0],new_angle[1],new_angle[2],*cost,orig_environment,*environment);
      printf("\n");
      }

   return(valid_transform);
   }

int valid_angle(angle,environment,b,t)
/*============================================================================*/
/* PURPOSE: LOOK UP THE SPECIFIED ANGLE IN THE LIST OF KNOWN
/*	ANGLES; IF IT EXISTS, RETURN THE ASSOCIATED PARAMETERS.
*/
   int angle[3];	/* ANGLE TO LOOK UP. */
   int* environment;	/* ANGLE ENVIRONMENT (4 MEMBERED RING, ETC). */
   double* b;		/* B PARAMETER OF MATCHED ANGLE, IF ANY. */
   double* t;		/* T PARAMETER OF MATCHED ANGLE, IF ANY. */
   {
   int i;		/* LOOP INDEX. */

   if (*environment==THREE_RING)
      {
      /* LOOK FOR THE ANGLE IN THE LIST OF KNOWN 3-MEMBERED RING ANGLES. */
      for (i=0; i<nangle_params_3ring; i++)
         {
         if ((angle_params_3ring[i].atom_types[0]==angle[0] &&
            angle_params_3ring[i].atom_types[1]==angle[1] &&
            angle_params_3ring[i].atom_types[2]==angle[2]) ||
            (angle_params_3ring[i].atom_types[0]==angle[2] &&
            angle_params_3ring[i].atom_types[1]==angle[1] &&
            angle_params_3ring[i].atom_types[2]==angle[0]))
            {
            /* FOUND THE SPECIFIED ANGLE, COPY PARAMETERS. */
            *b=angle_params_3ring[i].b;
            *t=angle_params_3ring[i].t;
            return(TRUE);
            }
         }
      }

   if (*environment==FOUR_RING)
      {
      /* LOOK FOR THE ANGLE IN THE LIST OF KNOWN 4-MEMBERED RING ANGLES. */
      for (i=0; i<nangle_params_4ring; i++)
         {
         if ((angle_params_4ring[i].atom_types[0]==angle[0] &&
            angle_params_4ring[i].atom_types[1]==angle[1] &&
            angle_params_4ring[i].atom_types[2]==angle[2]) ||
            (angle_params_4ring[i].atom_types[0]==angle[2] &&
            angle_params_4ring[i].atom_types[1]==angle[1] &&
            angle_params_4ring[i].atom_types[2]==angle[0]))
            {
            /* FOUND THE SPECIFIED ANGLE, COPY PARAMETERS. */
            *b=angle_params_4ring[i].b;
            *t=angle_params_4ring[i].t;
            return(TRUE);
            }
         }
      }

   if (*environment==NORMAL || *environment==THREE_RING || *environment==FOUR_RING)
      {
      /* LOOK FOR THE ANGLE IN THE LIST OF KNOWN NORMAL ANGLES. */
      for (i=0; i<nangle_params_normal; i++)
         {
         if ((angle_params_normal[i].atom_types[0]==angle[0] &&
            angle_params_normal[i].atom_types[1]==angle[1] &&
            angle_params_normal[i].atom_types[2]==angle[2]) ||
            (angle_params_normal[i].atom_types[0]==angle[2] &&
            angle_params_normal[i].atom_types[1]==angle[1] &&
            angle_params_normal[i].atom_types[2]==angle[0]))
            {
            /* FOUND THE SPECIFIED ANGLE, COPY PARAMETERS. */
            *b=angle_params_normal[i].b;
            *t=angle_params_normal[i].t;
            *environment=NORMAL;	/* MAY HAVE CHANGED ENVIRONMENT. */
            return(TRUE);
            }
         }
      }

   /* NO MATCH FOUND IN THE SPECIFIED ENVIRONMENT. */
   return(FALSE);
   }

int process_outofplane_errors()
/*============================================================================*/
/* PURPOSE: GENERATE MISSING OUT OF PLANE PARAMETERS BASED ON SIMILARITIES
/*	WITH EXISTING ATOM TYPE COMBINATIONS, OR ASK USER IF INSUFFICIENT
/*	SIMILARITIES EXIST.
*/
   {
   int outofplane[3];		/* ATOM TYPES IN OUT OF PLANE ANGLE. */
   double default_b;		/* BENDING CONSTANT PARAMETER. */
   int match;			/* TYPE OF MATCH FOUND FOR THIS PARAMETER. */
   char prompt[80];		/* PROMPT STRING TO USER. */
   int i;			/* LOOP INDEX. */
   int status;			/* USER IO STATUS. */

   for (i=0; i<noutofplane_errors; i++)
      {
      /* DEFINE THE ORIGINAL OUT OF PLANE, AND LOOK FOR A SIMILAR OUT OF PLANE. */
      outofplane[0]=outofplane_errors[i].atom_types[0];
      outofplane[1]=outofplane_errors[i].atom_types[1];
      outofplane[2]=outofplane_errors[i].atom_types[2];
      match=similar_outofplane_exists(outofplane,&default_b);

      /* PARAMETER VALUE WILL BE SAVED ONLY IF *NOT* AN EXACT MATCH. */
      outofplane_errors[i].save_param=(match==SIMILAR_MATCH || match==NO_MATCH);

      if (auto_take || match==EXACT_MATCH)
         {
         /* AUTOMATICALLY TAKE DEFAULT. */
         outofplane_errors[i].b=default_b;
         }
      else
         {
         /* ALLOW USER TO EDIT ONLY NON IDENTICALLY MATCHED PARAMETERS. */
         /* TELL THE USER WHICH OUT OF PLANE ANGLE WE ARE WORKING WITH. */
         printf("\nOut of plane %d-%d-%d",
            outofplane_errors[i].atom_types[0],
            outofplane_errors[i].atom_types[1],
            outofplane_errors[i].atom_types[2]);
         if (match==SIMILAR_MATCH)
            printf(":   (default from similar out of plane %d-%d-%d)\n",
               outofplane[0],outofplane[1],outofplane[2]);
         else
            printf(":   (no similar out of plane exists)\n");

         /* RUN THE DEFAULT PARAMETER BY THE USER FOR APPROVAL OR CHANGE. */
         sprintf(prompt,"Enter bending constant [%.3f]",default_b);
         status=gdouble(prompt,NULL_OK,RANGE_CHECK,0.0,1.0,&outofplane_errors[i].b);
         if (status==ABORT_ENTRY) return(FALSE);
         if (status==NULL_ENTRY) outofplane_errors[i].b=default_b;
         }
      }

   return(TRUE);
   }

int similar_outofplane_exists(outofplane,default_b)
/*============================================================================*/
/* PURPOSE: DETERMINE IF A SIMILAR OUT OF PLANE EXISTS.  IF SO, THE
/*	OUT OF PLANE CONSTANTS ARE RETRIEVED FROM THE SIMILAR DATA SET; IF NOT,
/*	THE OUT OF PLANE CONSTANTS ARE SET TO DEFAULT VALUES.  NOTE THAT THIS
/*	IS A SPECIAL CASE IN THAT ONLY ATOMS 1 AND 2 MAY BE TRANSFORMED (INDEED,
/*	ATOM 0 IS *ALWAYS* ZERO (AS A FLAG IN MM2 I BELIEVE)).
*/
   int outofplane[3];		/* ATOM TYPES IN OUT OF PLANE ANGLE. */
   double* default_b;		/* BENDING CONSTANT PARAMETER. */
   {
   int new_outofplane[3][3];	/* OUT OF PLANE WITH NTH ATOM TRANSFORMED. */
   int valid_transform[3];	/* FLAG; TRANSFORM OF NTH ATOM IS DEFINED. */
   int cost[3];			/* COST OF ACHIEVING NEW_OUTOFPLANE[N]. */
   int most_similar;		/* WHICH TRANSFORM IS MOST SIMILAR TO ORIG. */
   int i;			/* LOOP INDEX. */
   int j;			/* LOOP INDEX. */

   /* SET UP THE LAST DITCH DEFAULT IN CASE NO SIMILAR OUT OF PLANE EXISTS. */
   *default_b=0.05;

   /* SINCE WE KNOW MORE PARAMETERS THAN MM2, CHECK FOR AN EXACT MATCH FIRST. */
   if (valid_outofplane(outofplane,default_b))
      {
      if (debug_on) printf("\n   Out of plane %d-%d-%d found exact match.\n",outofplane[0],outofplane[1],outofplane[2]);
      return(EXACT_MATCH);
      }

   /* NO EXACT MATCH, SO TRY TRANSFORMING ONE ATOM IN EACH OF TWO COPIES OF */
   /* OUTOFPLANE, THEN PICK THE NEW OUTOFPLANE MOST SIMILAR TO THE ORIGINAL OUTOFPLANE. */
   if (debug_on) printf("\n   Out of plane %d-%d-%d transformations:",outofplane[0],outofplane[1],outofplane[2]);
   for (i=0; i<3; i++)
      {
      for (j=0; j<3; j++)
         {
         new_outofplane[i][j]=outofplane[j];
         }
      valid_transform[i]=find_valid_outofplane(outofplane,i,new_outofplane[i],&cost[i]);
      }

   /* IF BOTH ATOMS SUCCESSFULLY TRANSFORMED, PICK THE BEST. */
   if (valid_transform[1] && valid_transform[2])
      {
      if (cost[1] < cost[2]) most_similar=1;
      else if (cost[2] < cost[1]) most_similar=2;
      else
         {
         /* BOTH EXIST AND COST THE SAME, PICK LOWEST ORIGINAL ATOM TYPE. */
         if (outofplane[1] <= outofplane[2]) most_similar=1;
         else most_similar=2;
         }
      }
   else if (valid_transform[1]) most_similar=1;	/* ONLY ATOM 1 TRANSFORMED. */
   else if (valid_transform[2]) most_similar=2;	/* ONLY ATOM 2 TRANSFORMED. */
   else return(NO_MATCH);	/* NEITHER ATOM TRANSFORMED! */

   /* IF WE DROPPED THROUGH TO HERE, WE HAVE A VALID SIMILAR OUTOFPLANE, SO */
   /* COPY THE SINGLE CHANGED ATOM TYPE INTO THE ORIGINAL OUTOFPLANE AND */
   /* CALL VALID_OUTOFPLANE SOLELY TO RETRIEVE THE DEFAULT PARAMETERS. */
   outofplane[most_similar]=new_outofplane[most_similar][most_similar];
   valid_outofplane(outofplane,default_b);
   return(SIMILAR_MATCH);
   }

int find_valid_outofplane(orig_outofplane,which_atom,new_outofplane,cost)
/*============================================================================*/
/* PURPOSE: REPEATEDLY TRANSFORM THE SELECTED ATOM OF THE ORIGINAL OUT OF
/*	PLANE UNTIL WE EITHER FIND A VALID OUT OF PLANE FOR WHICH PARAMETERS
/*	EXIST OR WE RUN OUT OF TRANSFORMATIONS.
*/
   int orig_outofplane[3];	/* ORIGINAL OUT OF PLANE ATOM TYPES. */
   int which_atom;		/* WHICH ATOM TO TRANSFORM IN OUT OF PLANE. */
   int new_outofplane[3];	/* NEW, TRANSFORMED OUT OF PLANE IF SUCCESSFUL. */
   int* cost;			/* COST OF TRANSFORMATION. */
   {
   int transform_level;		/* WHICH TRANSFORMATION LEVEL TO USE. */
   int valid_transform;		/* FLAG; WHETHER A VALID TRANSFORM WAS FOUND. */
   double b;			/* B PARAMETER OF NEW OUT OF PLANE (JUNK HERE). */

   /* KEEP TRANSFORMING THE SELECTED ATOM UNTIL WE EITHER FIND A VALID */
   /* BOND FOR WHICH PARAMETERS EXIST OR WE RUN OUT OF TRANSFORMATIONS. */
   transform_level=0;	/* ALWAYS START AT THE BEGINNING. */
   while ((valid_transform=transform(orig_outofplane[which_atom],&transform_level,&new_outofplane[which_atom],cost)) &&
      !valid_outofplane(new_outofplane,&b))
      {
      if (debug_on)
         {
         printf("\n   Transforming atom #%d yields %2d-%2d-%2d: no parameters.",
            which_atom+1,new_outofplane[0],new_outofplane[1],new_outofplane[2]);
         }
      }

   if (debug_on)
      {
      if (valid_transform) printf("\n   Transforming atom #%d yields %2d-%2d-%2d: parameters exist (cost=%d).",
         which_atom+1,new_outofplane[0],new_outofplane[1],new_outofplane[2],*cost);
      printf("\n");
      }

   return(valid_transform);
   }

int valid_outofplane(outofplane,b)
/*============================================================================*/
/* PURPOSE: LOOK UP THE SPECIFIED OUT OF PLANE IN THE LIST OF KNOWN
/*	OUT OF PLANES; IF IT EXISTS, RETURN THE ASSOCIATED PARAMETERS.
*/
   int outofplane[3];	/* OUT OF PLANE TO LOOK UP. */
   double* b;		/* B PARAMETER OF MATCHED OUT OF PLANE, IF ANY. */
   {
   int i;		/* LOOP INDEX. */

   /* LOOK FOR THE SPECIFIED OUT OF PLANE IN THE LIST OF KNOWN OUT OF PLANES. */
   for (i=0; i<noutofplane_params; i++)
      {
      /* NOTE THAT OUT OF PLANE TYPES MUST MATCH EXACT ORDER. */
      if (outofplane_params[i].atom_types[0]==outofplane[0] &&
         outofplane_params[i].atom_types[1]==outofplane[1] &&
         outofplane_params[i].atom_types[2]==outofplane[2])
         {
         /* FOUND THE SPECIFIED OUT OF PLANE, COPY PARAMETERS. */
         *b=outofplane_params[i].b;
         return(TRUE);
         }
      }

   /* LOOKED THROUGH ALL KNOWN OUT OF PLANES, NO MATCH FOUND. */
   return(FALSE);
   }

int transform(original_atom_type,transform_level,new_atom_type,cost)
/*============================================================================*/
/* PURPOSE: TRANSFORM THE SPECIFIED ATOM BY THE TRANSFORM LEVEL
/*	ACCORDING TO THE LIST OF VALID TRANSFORMATIONS.  BOTH THE NEW ATOM
/*	TYPE AND THE COST OF THE COST OF THE TRANSFORMATION IS RETURNED.
/*	THE FUNCTION VALUE IS THE SUCCESS OF THE TRANSFORMATION (WHICH CAN ONLY
/*	FAIL IF THERE ARE NO MORE VALID TRANSFORMATIONS FOR THIS ATOM).  ON THE
/*	VERY FIRST TRANSFORMATION, TRANSFORM_LEVEL MUST BE 0 (TRANSFORM LEVEL
/*	IS INCREMENTED AT THE END OF THE ROUTINE TO BE READY FOR THE NEXT CALL).
*/
   int original_atom_type;	/* ORIGINAL ATOM TYPE BEFORE ANY TRANSFORMS. */
   int* transform_level;	/* NEW TRANSFORM LEVEL TO USE. */
   int* new_atom_type;		/* TRANSFORMED ATOM TYPE. */
   int* cost;			/* COST OF THE TRANSFORMATION. */
   {
   /* DEFINE THE TRANSFORMATIONS FOR SIMILAR ATOM TYPES.  FOR EACH ATOM */
   /* TYPE THERE IS A LIST OF UP TO MAX_TRANSFORMS TRANSFORMATIONS.  EACH */
   /* TRANSFORMATION DEFINES A PAIRED NEW ATOM TYPE AND A COST. */
   static int similar[MAX_ATOM_TYPES][MAX_TRANSFORMS][2] =
      {
      { {5,1}                  },	/* ATOM TYPE 1 TRANSFORMATIONS. */
      { {1,5},  {5,10}         },	/* ATOM TYPE 2 TRANSFORMATIONS. */
      { {2,1}                  },	/* ATOM TYPE 3 TRANSFORMATIONS. */
      { {2,2},  {1,10}, {5,11} },	/* ATOM TYPE 4 TRANSFORMATIONS. */
      { {0,0}                  },	/* ATOM TYPE 5 TRANSFORMATIONS. */
      { {8,1},  {1,5}          },	/* ATOM TYPE 6 TRANSFORMATIONS. */
      { {2,5}                  },	/* ATOM TYPE 7 TRANSFORMATIONS. */
      { {6,1},  {1,5}          },	/* ATOM TYPE 6 TRANSFORMATIONS. */
      { {8,1},  {6,2}          },	/* ATOM TYPE 9 TRANSFORMATIONS. */
      { {4,1},  {1,10}, {5,11} },	/* ATOM TYPE 10 TRANSFORMATIONS. */
      { {5,1},  {12,1}         },	/* ATOM TYPE 11 TRANSFORMATIONS. */
      { {5,1}                  },	/* ATOM TYPE 12 TRANSFORMATIONS. */
      { {12,1}, {5,1}          },	/* ATOM TYPE 13 TRANSFORMATIONS. */
      { {13,1}, {12,1}, {5,1}  },	/* ATOM TYPE 14 TRANSFORMATIONS. */
      { {6,1},  {1,5}          },	/* ATOM TYPE 15 TRANSFORMATIONS. */
      { {15,1}, {6,2},  {1,5}  },	/* ATOM TYPE 16 TRANSFORMATIONS. */
      { {3,1},  {2,2}          },	/* ATOM TYPE 17 TRANSFORMATIONS. */
      { {3,1},  {2,2}          },	/* ATOM TYPE 18 TRANSFORMATIONS. */
      { {1,1}                  },	/* ATOM TYPE 19 TRANSFORMATIONS. */
      { {0,0}                  },	/* ATOM TYPE 20 TRANSFORMATIONS. */
      { {23,1}, {5,1}          },	/* ATOM TYPE 21 TRANSFORMATIONS. */
      { {1,1}                  },	/* ATOM TYPE 22 TRANSFORMATIONS. */
      { {21,1}, {5,1}          },	/* ATOM TYPE 23 TRANSFORMATIONS. */
      { {21,1}, {5,1}          },	/* ATOM TYPE 24 TRANSFORMATIONS. */
      { {0,0}                  },	/* ATOM TYPE 25 TRANSFORMATIONS. */
      { {0,0}                  },	/* ATOM TYPE 26 TRANSFORMATIONS. */
      { {0,0}                  },	/* ATOM TYPE 27 TRANSFORMATIONS. */
      { {23,1}, {21,1}, {5,1}  },	/* ATOM TYPE 28 TRANSFORMATIONS. */
      { {0,0}                  },	/* ATOM TYPE 29 TRANSFORMATIONS. */
      { {0,0}                  },	/* ATOM TYPE 30 TRANSFORMATIONS. */
      { {0,0}                  },	/* ATOM TYPE 31 TRANSFORMATIONS. */
      { {0,0}                  },	/* ATOM TYPE 32 TRANSFORMATIONS. */
      { {0,0}                  },	/* ATOM TYPE 33 TRANSFORMATIONS. */
      { {0,0}                  },	/* ATOM TYPE 34 TRANSFORMATIONS. */
      { {0,0}                  },	/* ATOM TYPE 35 TRANSFORMATIONS. */
      { {5,1}                  },	/* ATOM TYPE 36 TRANSFORMATIONS. */
      { {40,1}, {2,1}          },	/* ATOM TYPE 37 TRANSFORMATIONS. */
      { {2,1}                  },	/* ATOM TYPE 38 TRANSFORMATIONS. */
      { {8,1},  {1,5}          },	/* ATOM TYPE 39 TRANSFORMATIONS. */
      { {37,1}, {2,1}          },	/* ATOM TYPE 40 TRANSFORMATIONS. */
      { {37,1}, {2,1}          },	/* ATOM TYPE 41 TRANSFORMATIONS. */
      { {41,1}, {2,2}          },	/* ATOM TYPE 42 TRANSFORMATIONS. */
      { {37,1}, {2,2}          },	/* ATOM TYPE 43 TRANSFORMATIONS. */
      { {5,0},  {21,1}         },	/* ATOM TYPE 44 TRANSFORMATIONS. */
      { {43,1}, {37,2}, {2,5}  },	/* ATOM TYPE 45 TRANSFORMATIONS. */
      { {3,2},  {2,5}          },	/* ATOM TYPE 46 TRANSFORMATIONS. */
      { {6,1}                  },	/* ATOM TYPE 47 TRANSFORMATIONS. */
      { {23,1}, {5,1}          },	/* ATOM TYPE 48 TRANSFORMATIONS. */
      { {6,1},  {22,5}, {1,6}  },	/* ATOM TYPE 49 TRANSFORMATIONS. */
      { {2,0}                  },	/* ATOM TYPE 50 TRANSFORMATIONS. */
      { {0,0}                  },	/* ATOM TYPE 51 TRANSFORMATIONS. */
      { {0,0}                  },	/* ATOM TYPE 52 TRANSFORMATIONS. */
      { {0,0}                  },	/* ATOM TYPE 53 TRANSFORMATIONS. */
      { {0,0}                  },	/* ATOM TYPE 54 TRANSFORMATIONS. */
      { {0,0}                  },	/* ATOM TYPE 55 TRANSFORMATIONS. */
      { {1,0}                  },	/* ATOM TYPE 56 TRANSFORMATIONS. */
      { {2,0},  {1,5}          },	/* ATOM TYPE 57 TRANSFORMATIONS. */
      { {3,0},  {2,1}          },	/* ATOM TYPE 58 TRANSFORMATIONS. */
      { {0,0}                  },	/* ATOM TYPE 59 TRANSFORMATIONS. */
      { {25,0}                 }	/* ATOM TYPE 60 TRANSFORMATIONS. */
      };

   /* CHECK THAT THIS TRANSFORMATION REQUEST IS VALID. */
   if (*transform_level>=MAX_TRANSFORMS || original_atom_type<1) return(FALSE);
   else
      {
      /* VALID TRANSFORMATION REQUESTED; RETURN NEW TYPE AND COST IF DEFINED. */
      *new_atom_type=similar[original_atom_type-1][*transform_level][ATOM_TYPE];
      if (*new_atom_type==0) return(FALSE);
      else
         {
         /* THIS IS A DEFINED TRANSFORMATION; GET THE COST AND INCREMENT THE LEVEL FOR NEXT TIME. */
         *cost=similar[original_atom_type-1][*transform_level][COST];
         (*transform_level)++;
         return(TRUE);
         }
      }
   }

void write_parameters(param_filename)
/*============================================================================*/
/* PURPOSE: WRITE A CORRECTLY FORMATTED MM2/MM3 PARAMETER FILE CONTAINING ALL
/*	THE NEW PARAMETERS.  LINE DESCRIPTIONS ARE VERBATIM FROM THE MM2(87)
/*	QCPE PROGRAM MANUAL.
*/
   char param_filename[];	/* NEW PARAMETER FILE NAME. */
   {
   FILE* param_file;		/* NEW PARAMETER FILE. */
   char env_flag;		/* ORIGINAL ENVIRONMENT FLAG (LINES 9A,9F). */
   int i;			/* LOOP INDEX. */

   param_file=fopen(param_filename,"w");

   /* WRITE LINE 9 (CHANGED CONSTANTS). */
   fprintf(param_file,"%5d%5d%2d%3d%5d%5d%5d%5d%5d\n",
      ntorsion_errors,nbond_errors,0,0,nangle_errors+noutofplane_errors,nbond_errors,0,0,0);

   /* WRITE LINES 9A (TORSIONAL CONSTANTS). */
   for (i=0; i<ntorsion_errors; i++)
      {
      if (torsion_errors[i].orig_env==FOUR_RING) env_flag='4';
      else env_flag=' ';

      fprintf(param_file,"%1c%4d%5d%5d%5d%10.5f%10.5f%10.5f\n",
         env_flag,
         torsion_errors[i].atom_types[0],torsion_errors[i].atom_types[1],
         torsion_errors[i].atom_types[2],torsion_errors[i].atom_types[3],
         torsion_errors[i].v1,torsion_errors[i].v2,torsion_errors[i].v3);
      }

   /* WRITE LINES 9B (STRETCHING AND ANOMERIC EFFECT CONSTANTS). */
   for (i=0; i<nbond_errors; i++)
      {
      fprintf(param_file,"%5d%5d%10.5f%10.5f%10.5f%10.5f%10.5f\n",
         bond_errors[i].atom_types[0],bond_errors[i].atom_types[1],
         bond_errors[i].s,bond_errors[i].l0,0.0,bond_errors[i].slps,bond_errors[i].slpt);
      }

   /* WRITE LINES 9C (BOND DIPOLE CONSTANTS BY ATOM TYPE). */
   for (i=0; i<nbond_errors; i++)
      {
      fprintf(param_file,"%5d%5d%10.5f\n",
         bond_errors[i].atom_types[0],bond_errors[i].atom_types[1],
         bond_errors[i].bmom);
      }

   /* WRITE LINES 9F (BENDING CONSTANTS); NOTE THAT ALL NORMAL ANGLES */
   /* (FIRST ATOM TYPE NOT 0) COME FIRST, THEN ALL OUT-OF-PLANE VALUES */
   /* (FIRST ATOM TYPE 0) WHICH ONLY HAVE BENDING CONSTANT. */
   for (i=0; i<nangle_errors; i++)
      {
      if (angle_errors[i].orig_env==THREE_RING) env_flag='3';
      else if (angle_errors[i].orig_env==FOUR_RING) env_flag='4';
      else env_flag=' ';

      fprintf(param_file,"%1c%4d%5d%5d     %10.5f%10.5f    0\n",
         env_flag,
         angle_errors[i].atom_types[0],angle_errors[i].atom_types[1],angle_errors[i].atom_types[2],
         angle_errors[i].b,angle_errors[i].t);
      }
   for (i=0; i<noutofplane_errors; i++)
      {
      fprintf(param_file," %4d%5d%5d     %10.5f\n",
         outofplane_errors[i].atom_types[0],outofplane_errors[i].atom_types[1],outofplane_errors[i].atom_types[2],
         outofplane_errors[i].b);
      }

   fclose(param_file);

   return;
   }