
#include "msi2lmp.h"
#include "Forcefield.h"

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

static int find_improper_body_data(char [][5],struct FrcFieldItem,int *);
static void rearrange_improper(int,int);
static int find_trigonal_body_data(char [][5],struct FrcFieldItem);
static int find_angleangle_data(char [][5],struct FrcFieldItem,int[]);
static int find_match(int, char [][5],struct FrcFieldItem,int *);
static int match_types(int,int,char [][5],char [][5],int *);
static double get_r0(int,int);
static double get_t0(int,int,int);
static int quo_cp();
static void get_equivs(int,char [][5],char[][5]);
static int find_equiv_type(char[]);

/**********************************************************************/
/*                                                                    */
/*  GetParameters is a long routine for searching the forcefield      */
/*  parameters (read in by ReadFrcFile) for parameters corresponding  */
/*  to the different internal coordinate types derived by MakeLists   */
/*                                                                    */
/**********************************************************************/

void GetParameters()
{
  int i,j,k,backwards,cp_type,rearrange;
  int kloc[3],multiplicity;
  char potential_types[4][5];
  char equiv_types[4][5];
  double rab,rbc,rcd,tabc,tbcd,tabd,tcbd;

  if (pflag > 1) fprintf(stderr," Trying Atom Equivalences if needed\n");

  /**********************************************************************/
  /*                                                                    */
  /*   Find masses  of atom types                                       */
  /*                                                                    */
  /**********************************************************************/

  for (i=0; i < no_atom_types; i++) {
    backwards = -1;
    strncpy(potential_types[0],atomtypes[i].potential,5);
    k = find_match(1,potential_types,ff_atomtypes,&backwards);
    if (k < 0) {
      printf(" Unable to find mass for %s\n",atomtypes[i].potential);
      condexit(10);
    } else {
      atomtypes[i].mass = ff_atomtypes.data[k].ff_param[0];
    }
  }

  /**********************************************************************/
  /*                                                                    */
  /*   Find VDW parameters for atom types                               */
  /*                                                                    */
  /**********************************************************************/

  for (i=0; i < no_atom_types; i++) {
    backwards = 0;
    for (j=0; j < 2; j++) atomtypes[i].params[j] = 0.0;
    strncpy(potential_types[0],atomtypes[i].potential,5);
    k = find_match(1,potential_types,ff_vdw,&backwards);
    if (k < 0) {
      get_equivs(1,potential_types,equiv_types);

      if (pflag > 2) printf(" Using equivalences for VDW %s -> %s\n",
                            potential_types[0],equiv_types[0]);

      k = find_match(1,equiv_types,ff_vdw,&backwards);
    }
    if (k < 0) {
      printf(" Unable to find vdw data for %s\n",atomtypes[i].potential);
      condexit(11);
    } else {
      if (ljtypeflag == 0) {
        if((ff_vdw.data[k].ff_param[0] != 0.0 ) &&
           (ff_vdw.data[k].ff_param[1] != 0.0)) {
          atomtypes[i].params[0] =
            (ff_vdw.data[k].ff_param[1]*
             ff_vdw.data[k].ff_param[1])/(4.0*ff_vdw.data[k].ff_param[0]);
          atomtypes[i].params[1] = pow((ff_vdw.data[k].ff_param[0]/
                                        ff_vdw.data[k].ff_param[1]),
                                       (1.0/6.0));
        }
      } else if (ljtypeflag == 1) {
        atomtypes[i].params[0] = ff_vdw.data[k].ff_param[1];
        atomtypes[i].params[1] = ff_vdw.data[k].ff_param[0];
      } else {
          printf(" Unknown LJ parameter type %d\n",ljtypeflag);
          exit(111);
      }
    }
  }

  if (pflag > 2) {
    printf("\n Atom Types, Masses and VDW Parameters\n");
    for (i=0; i < no_atom_types; i++) {
      printf(" %3s %8.4f %8.4f %8.4f\n",
             atomtypes[i].potential,atomtypes[i].mass, atomtypes[i].params[0],atomtypes[i].params[1]);
    }
  }

  /**********************************************************************/
  /*                                                                    */
  /*   Find parameters for bond types                                   */
  /*                                                                    */
  /**********************************************************************/

  for (i=0; i < no_bond_types; i++) {
    backwards = 0;
    for (j=0; j < 4; j++) bondtypes[i].params[j] = 0.0;
    for (j=0; j < 2; j++)
      strncpy(potential_types[j],
              atomtypes[bondtypes[i].types[j]].potential,5);
    k = find_match(2,potential_types,ff_bond,&backwards);
    if (k < 0) {
      get_equivs(2,potential_types,equiv_types);

      if (pflag > 2) {
        printf(" Using equivalences for bond %s %s -> %s %s\n",
               potential_types[0],potential_types[1],
               equiv_types[0],equiv_types[1]);
      }
      k = find_match(2,equiv_types,ff_bond,&backwards);
    }
    if (k < 0) {
      printf(" Unable to find bond data for %s %s\n",
             potential_types[0],potential_types[1]);
      condexit(12);
    } else {
      if (forcefield & (FF_TYPE_CLASS1|FF_TYPE_OPLSAA)) {
        bondtypes[i].params[0] = ff_bond.data[k].ff_param[1];
        bondtypes[i].params[1] = ff_bond.data[k].ff_param[0];
      }

      if (forcefield & FF_TYPE_CLASS2) {
        for (j=0; j < 4; j++)
          bondtypes[i].params[j] = ff_bond.data[k].ff_param[j];
      }
    }
  }

  if (pflag > 2) {
    printf("\n Bond Types and  Parameters\n");
    for (i=0; i < no_bond_types; i++) {
      for (j=0; j < 2; j++)
        printf(" %-3s",atomtypes[bondtypes[i].types[j]].potential);
      for (j=0; j < 4; j++)
        printf(" %8.4f",bondtypes[i].params[j]);
      printf("\n");
    }
  }


  /**********************************************************************/
  /*                                                                    */
  /*   Find parameters for angle types including bondbond,              */
  /*   and bondangle parameters if Class II                             */
  /*                                                                    */
  /*   Each of the cross terms are searched separately even though      */
  /*   they share a given angle type. This allows parameters to be      */
  /*   in different order in the forcefield for each cross term or      */
  /*   maybe not even there.                                            */
  /*                                                                    */
  /**********************************************************************/
  for (i=0; i < no_angle_types; i++) {
    backwards = 0;
    for (j=0; j < 4; j++) angletypes[i].params[j] = 0.0;
    for (j=0; j < 3; j++)
      strncpy(potential_types[j],atomtypes[angletypes[i].types[j]].potential,5);
    k = find_match(3,potential_types,ff_ang,&backwards);
    if (k < 0) {
      get_equivs(3,potential_types,equiv_types);
      if (pflag > 2) {
        printf(" Using equivalences for angle %s %s %s -> %s %s %s\n",
               potential_types[0],potential_types[1],
               potential_types[2],
               equiv_types[0],equiv_types[1],
               equiv_types[2]);
      }
      k = find_match(3,equiv_types,ff_ang,&backwards);
    }
    if (k < 0) {
      printf(" Unable to find angle data for %s %s %s\n",
             potential_types[0],potential_types[1],potential_types[2]);
      condexit(13);
    } else {
      if (forcefield & (FF_TYPE_CLASS1|FF_TYPE_OPLSAA)) {
        angletypes[i].params[0] = ff_ang.data[k].ff_param[1];
        angletypes[i].params[1] = ff_ang.data[k].ff_param[0];
      }

      if (forcefield & FF_TYPE_CLASS2) {
        for (j=0; j < 4; j++)
          angletypes[i].params[j] = ff_ang.data[k].ff_param[j];
      }
    }
    if (forcefield & FF_TYPE_CLASS2) {
      get_equivs(3,potential_types,equiv_types);
      if (pflag > 2) {
        printf(" Using equivalences for 3 body cross terms %s %s %s -> %s %s %s\n",
               potential_types[0],potential_types[1],potential_types[2],
               equiv_types[0],equiv_types[1],equiv_types[2]);
      }
      for (j=0; j < 3; j++) angletypes[i].bondbond_cross_term[j] = 0.0;
      for (j=0; j < 4; j++) angletypes[i].bondangle_cross_term[j] = 0.0;

      rab = get_r0(angletypes[i].types[0],angletypes[i].types[1]);
      rbc = get_r0(angletypes[i].types[1],angletypes[i].types[2]);

      angletypes[i].bondbond_cross_term[1] = rab;
      angletypes[i].bondbond_cross_term[2] = rbc;
      angletypes[i].bondangle_cross_term[2] = rab;
      angletypes[i].bondangle_cross_term[3] = rbc;

      k = find_match(3,potential_types,ff_bonbon,&backwards);
      if (k < 0) {
        k = find_match(3,equiv_types,ff_bonbon,&backwards);
      }
      if (k < 0) {
        printf(" Unable to find bondbond data for %s %s %s\n",
               potential_types[0],potential_types[1],potential_types[2]);
        condexit(14);
      } else {
        angletypes[i].bondbond_cross_term[0] = ff_bonbon.data[k].ff_param[0];
      }
      k = find_match(3,potential_types,ff_bonang,&backwards);
      if (k < 0) {
        k = find_match(3,equiv_types,ff_bonang,&backwards);
      }
      if (k < 0) {
        printf(" Unable to find bondangle data for %s %s %s\n",
               potential_types[0],potential_types[1],potential_types[2]);
        condexit(15);
      } else {
        if (backwards) {
          angletypes[i].bondangle_cross_term[0] = ff_bonang.data[k].ff_param[1];
          angletypes[i].bondangle_cross_term[1] = ff_bonang.data[k].ff_param[0];
        } else {
          angletypes[i].bondangle_cross_term[0] = ff_bonang.data[k].ff_param[0];
          angletypes[i].bondangle_cross_term[1] = ff_bonang.data[k].ff_param[1];
        }
      }
    }
  }

  if (pflag > 2) {
    printf("\n Angle Types and Parameters\n");
    for (i=0; i < no_angle_types; i++) {
      for (j=0; j < 3; j++)
        printf(" %-3s", atomtypes[angletypes[i].types[j]].potential);
      for (j=0; j < 4; j++) printf(" %8.4f",angletypes[i].params[j]);
      printf("\n");
    }

    if (forcefield & FF_TYPE_CLASS2) {
      printf("\n BondBond Types and  Parameters\n");
      for (i=0; i < no_angle_types; i++) {
        for (j=0; j < 3; j++)
          printf(" %-3s",atomtypes[angletypes[i].types[j]].potential);
        for (j=0; j < 3; j++)
          printf(" %8.4f",angletypes[i].bondbond_cross_term[j]);
        printf("\n");
      }
      printf("\n BondAngle Types and  Parameters\n");
      for (i=0; i < no_angle_types; i++) {
        for (j=0; j < 3; j++)
          printf(" %-3s",atomtypes[angletypes[i].types[j]].potential);
        for (j=0; j < 4; j++)
          printf(" %8.4f",angletypes[i].bondangle_cross_term[j]);
        printf("\n");
      }
    }
  }

  /**********************************************************************/
  /*                                                                    */
  /*   Find parameters for dihedral types including endbonddihedral,    */
  /*   midbonddihedral, angledihedral, angleangledihedral and           */
  /*   bondbond13 parameters if Class II                                */
  /*                                                                    */
  /*   Each of the cross terms are searched separately even though      */
  /*   they share a given dihedral type. This allows parameters to be   */
  /*   in different order in the forcefield for each cross term or      */
  /*   maybe not even there.                                            */
  /*                                                                    */
  /**********************************************************************/

  for (i=0; i < no_dihedral_types; i++) {
    for (j=0; j < 6; j++)
      dihedraltypes[i].params[j] = 0.0;
    for (j=0; j < 4; j++)
      strncpy(potential_types[j],
              atomtypes[dihedraltypes[i].types[j]].potential,5);
    backwards = 0;
    k = find_match(4,potential_types,ff_tor,&backwards);

    if (k < 0) {
      get_equivs(4,potential_types,equiv_types);

      if (pflag > 2) {
        printf(" Using equivalences for dihedral %s %s %s %s -> %s %s %s %s\n",
               potential_types[0],potential_types[1],
               potential_types[2],potential_types[3],
               equiv_types[0],equiv_types[1],
               equiv_types[2],equiv_types[3]);
      }
      k = find_match(4,equiv_types,ff_tor,&backwards);
    }
    if (k < 0) {
      printf(" Unable to find torsion data for %s %s %s %s\n",
             potential_types[0],
             potential_types[1],
             potential_types[2],
             potential_types[3]);
      condexit(16);
    } else {
      if (forcefield & FF_TYPE_CLASS1) {
        multiplicity = 1;
        if (ff_tor.data[k].ff_types[0][0] == '*')
          multiplicity =
            atomtypes[dihedraltypes[i].types[1]].no_connect-1;
        if (ff_tor.data[k].ff_types[3][0] == '*')
          multiplicity *=
            atomtypes[dihedraltypes[i].types[2]].no_connect-1;

        dihedraltypes[i].params[0] = ff_tor.data[k].ff_param[0]/(double) multiplicity;
        if (ff_tor.data[k].ff_param[2] == 0.0)
          dihedraltypes[i].params[1] = 1.0;
        else if (ff_tor.data[k].ff_param[2] == 180.0)
          dihedraltypes[i].params[1] = -1.0;
        else {
          printf(" Non planar phi0 for %s %s %s %s\n",
                 potential_types[0],potential_types[1],
                 potential_types[2],potential_types[3]);
          dihedraltypes[i].params[1] = 0.0;
        }
        dihedraltypes[i].params[2] = ff_tor.data[k].ff_param[1];
      }
      if (forcefield & FF_TYPE_OPLSAA) {
        for (j=0; j < 4; j++)
          dihedraltypes[i].params[j] = ff_tor.data[k].ff_param[j];
      }
      if (forcefield & FF_TYPE_CLASS2) {
        for (j=0; j < 6; j++)
          dihedraltypes[i].params[j] = ff_tor.data[k].ff_param[j];
      }
    }

    if (forcefield & FF_TYPE_CLASS2) {
      get_equivs(4,potential_types,equiv_types);
      if (pflag > 2) {
        printf(" Using equivalences for linear 4 body cross terms  %s %s %s %s -> %s %s %s %s\n",
               potential_types[0],potential_types[1],
               potential_types[2],potential_types[3],
               equiv_types[0],equiv_types[1],
               equiv_types[2],equiv_types[3]);
      }

      for (j=0; j < 8; j++)
        dihedraltypes[i].endbonddihedral_cross_term[j] = 0.0;
      for (j=0; j < 4; j++)
        dihedraltypes[i].midbonddihedral_cross_term[j] = 0.0;
      for (j=0; j < 8; j++)
        dihedraltypes[i].angledihedral_cross_term[j] = 0.0;
      for (j=0; j < 3; j++)
        dihedraltypes[i].angleangledihedral_cross_term[j] = 0.0;
      for (j=0; j < 3; j++)
        dihedraltypes[i].bond13_cross_term[j] = 0.0;

      rab = get_r0(dihedraltypes[i].types[0],dihedraltypes[i].types[1]);
      rbc = get_r0(dihedraltypes[i].types[1],dihedraltypes[i].types[2]);
      rcd = get_r0(dihedraltypes[i].types[2],dihedraltypes[i].types[3]);
      tabc = get_t0(dihedraltypes[i].types[0],
                    dihedraltypes[i].types[1],
                    dihedraltypes[i].types[2]);

      tbcd = get_t0(dihedraltypes[i].types[1],
                    dihedraltypes[i].types[2],
                    dihedraltypes[i].types[3]);

      dihedraltypes[i].endbonddihedral_cross_term[6] = rab;
      dihedraltypes[i].endbonddihedral_cross_term[7] = rcd;
      dihedraltypes[i].midbonddihedral_cross_term[3] = rbc;
      dihedraltypes[i].angledihedral_cross_term[6] = tabc;
      dihedraltypes[i].angledihedral_cross_term[7] = tbcd;
      dihedraltypes[i].angleangledihedral_cross_term[1] = tabc;
      dihedraltypes[i].angleangledihedral_cross_term[2] = tbcd;
      dihedraltypes[i].bond13_cross_term[1] = rab;
      dihedraltypes[i].bond13_cross_term[2] = rcd;

      backwards = 0;
      k = find_match(4,potential_types,ff_endbontor,&backwards);
      if (k < 0) {
        k = find_match(4,equiv_types,ff_endbontor,&backwards);
      }
      if (k < 0) {
        printf(" Unable to find endbonddihedral data for %s %s %s %s\n",
               potential_types[0],potential_types[1],
               potential_types[2],potential_types[3]);
        condexit(17);
      } else {
        if (backwards) {
          dihedraltypes[i].endbonddihedral_cross_term[0] =
            ff_endbontor.data[k].ff_param[3];
          dihedraltypes[i].endbonddihedral_cross_term[1] =
            ff_endbontor.data[k].ff_param[4];
          dihedraltypes[i].endbonddihedral_cross_term[2] =
            ff_endbontor.data[k].ff_param[5];
          dihedraltypes[i].endbonddihedral_cross_term[3] =
            ff_endbontor.data[k].ff_param[0];
          dihedraltypes[i].endbonddihedral_cross_term[4] =
            ff_endbontor.data[k].ff_param[1];
          dihedraltypes[i].endbonddihedral_cross_term[5] =
            ff_endbontor.data[k].ff_param[2];
        } else {
          dihedraltypes[i].endbonddihedral_cross_term[0] =
            ff_endbontor.data[k].ff_param[0];
          dihedraltypes[i].endbonddihedral_cross_term[1] =
            ff_endbontor.data[k].ff_param[1];
          dihedraltypes[i].endbonddihedral_cross_term[2] =
            ff_endbontor.data[k].ff_param[2];
          dihedraltypes[i].endbonddihedral_cross_term[3] =
            ff_endbontor.data[k].ff_param[3];
          dihedraltypes[i].endbonddihedral_cross_term[4] =
            ff_endbontor.data[k].ff_param[4];
          dihedraltypes[i].endbonddihedral_cross_term[5] =
            ff_endbontor.data[k].ff_param[5];
        }
      }
      backwards = 0;
      k = find_match(4,potential_types,ff_midbontor,&backwards);
      if (k < 0) {
        k = find_match(4,equiv_types,ff_midbontor,&backwards);
      }
      if (k < 0) {
        printf(" Unable to find midbonddihedral data for %s %s %s %s\n",
               potential_types[0],potential_types[1],
               potential_types[2],potential_types[3]);
        condexit(18);
      } else {
        dihedraltypes[i].midbonddihedral_cross_term[0] =
          ff_midbontor.data[k].ff_param[0];
        dihedraltypes[i].midbonddihedral_cross_term[1] =
          ff_midbontor.data[k].ff_param[1];
        dihedraltypes[i].midbonddihedral_cross_term[2] =
          ff_midbontor.data[k].ff_param[2];
      }

      backwards = 0;
      k = find_match(4,potential_types,ff_angtor,&backwards);
      if (k < 0) {
        k = find_match(4,equiv_types,ff_angtor,&backwards);
      }
      if (k < 0) {
        printf(" Unable to find angledihedral data for %s %s %s %s\n",
               potential_types[0],potential_types[1],
               potential_types[2],potential_types[3]);
        condexit(19);
      } else {
        if (backwards) {
          dihedraltypes[i].angledihedral_cross_term[0] =
            ff_angtor.data[k].ff_param[3];
          dihedraltypes[i].angledihedral_cross_term[1] =
            ff_angtor.data[k].ff_param[4];
          dihedraltypes[i].angledihedral_cross_term[2] =
            ff_angtor.data[k].ff_param[5];
          dihedraltypes[i].angledihedral_cross_term[3] =
            ff_angtor.data[k].ff_param[0];
          dihedraltypes[i].angledihedral_cross_term[4] =
            ff_angtor.data[k].ff_param[1];
          dihedraltypes[i].angledihedral_cross_term[5] =
            ff_angtor.data[k].ff_param[2];
        } else {
          dihedraltypes[i].angledihedral_cross_term[0] =
            ff_angtor.data[k].ff_param[0];
          dihedraltypes[i].angledihedral_cross_term[1] =
            ff_angtor.data[k].ff_param[1];
          dihedraltypes[i].angledihedral_cross_term[2] =
            ff_angtor.data[k].ff_param[2];
          dihedraltypes[i].angledihedral_cross_term[3] =
            ff_angtor.data[k].ff_param[3];
          dihedraltypes[i].angledihedral_cross_term[4] =
            ff_angtor.data[k].ff_param[4];
          dihedraltypes[i].angledihedral_cross_term[5] =
            ff_angtor.data[k].ff_param[5];
        }
      }
      backwards = 0;
      k = find_match(4,potential_types,ff_angangtor,&backwards);
      if (k < 0) {
        k = find_match(4,equiv_types,ff_angangtor,&backwards);
      }
      if (k < 0) {
        printf(" Unable to find angleangledihedral data for %s %s %s %s\n",
               potential_types[0],potential_types[1],
               potential_types[2],potential_types[3]);
        condexit(20);
      } else {
        dihedraltypes[i].angleangledihedral_cross_term[0] =
          ff_angangtor.data[k].ff_param[0];
      }
      cp_type = quo_cp();
      if ((cp_type >= 0) &&
          ((dihedraltypes[i].types[0] == cp_type) ||
           (dihedraltypes[i].types[1] == cp_type) ||
           (dihedraltypes[i].types[2] == cp_type) ||
           (dihedraltypes[i].types[3] == cp_type)   )) {
        backwards = 0;
        k = find_match(4,potential_types,ff_bonbon13,&backwards);
        if (k < 0) {
          k = find_match(4,equiv_types,ff_bonbon13,&backwards);
        }
        if (k < 0) {
          printf(" Unable to find bond13 data for %s %s %s %s\n",
                 potential_types[0],potential_types[1],
                 potential_types[2],potential_types[3]);
          condexit(21);
        } else {
          dihedraltypes[i].bond13_cross_term[0] =
            ff_bonbon13.data[k].ff_param[0];
        }
      }
    }
  }

  if (pflag > 2) {
    printf("\n Dihedral Types and  Parameters\n");
    for (i=0; i < no_dihedral_types; i++) {
      for (j=0; j < 4; j++)
        printf(" %-3s",atomtypes[dihedraltypes[i].types[j]].potential);
      for (j=0; j < 6; j++)
        printf(" %8.4f",dihedraltypes[i].params[j]);
      printf("\n");
    }

    if (forcefield & FF_TYPE_CLASS2) {

      printf("\n EndBondDihedral Types and Parameters\n");
      for (i=0; i < no_dihedral_types; i++) {
        for (j=0; j < 4; j++)
          printf(" %-3s",atomtypes[dihedraltypes[i].types[j]].potential);
        for (j=0; j < 8; j++)
          printf(" %8.4f",dihedraltypes[i].endbonddihedral_cross_term[j]);
        printf("\n");
      }
      printf("\n MidBondDihedral Types and Parameters\n");
      for (i=0; i < no_dihedral_types; i++) {
        for (j=0; j < 4; j++)
          printf(" %-3s",atomtypes[dihedraltypes[i].types[j]].potential);
        for (j=0; j < 4; j++)
          printf(" %8.4f",dihedraltypes[i].midbonddihedral_cross_term[j]);
        printf("\n");
      }

      printf("\n AngleDihedral Types and Parameters\n");
      for (i=0; i < no_dihedral_types; i++) {
        for (j=0; j < 4; j++)
          printf(" %-3s",atomtypes[dihedraltypes[i].types[j]].potential);
        for (j=0; j < 8; j++)
          printf(" %8.4f",dihedraltypes[i].angledihedral_cross_term[j]);
        printf("\n");
      }

      printf("\n AngleAngleDihedral Types and Parameters\n");
      for (i=0; i < no_dihedral_types; i++) {
        for (j=0; j < 4; j++)
          printf(" %-3s",atomtypes[dihedraltypes[i].types[j]].potential);
        for (j=0; j < 3; j++)
          printf(" %8.4f",dihedraltypes[i].angleangledihedral_cross_term[j]);
        printf("\n");
      }

      printf("\n Bond13 Types and  Parameters\n");

      for (i=0; i < no_dihedral_types; i++) {
        for (j=0; j < 4; j++)
          printf(" %-3s",atomtypes[dihedraltypes[i].types[j]].potential);
        for (j=0; j < 3; j++)
          printf(" %8.4f",dihedraltypes[i].bond13_cross_term[j]);
        printf("\n");
      }
    }
  }


  /**********************************************************************/
  /*                                                                    */
  /*   Find parameters for oop types                                    */
  /*                                                                    */
  /*   This is the most complicated of all the types because the        */
  /*   the class I oop is actually an improper torsion and does         */
  /*   not have the permutation symmetry of a well defined oop          */
  /*   The net result is that if one does not find the current          */
  /*   atom type ordering in the forcefield file then one must try each */
  /*   of the next permutations (6 in total) and when a match is found  */
  /*   the program must go back and rearrange the oop type AND the atom */
  /*   ordering in the oop lists for those with the current type        */
  /*                                                                    */
  /*   The Class II oop types are easier but also tedious since the     */
  /*   program has to try all permutations of the a c and d atom        */
  /*   types to find a match. A special routine is used to do this.     */
  /*                                                                    */
  /*   Fortunately, there are typically few oop types                   */
  /*                                                                    */
  /**********************************************************************/

  if (forcefield & FF_TYPE_CLASS1) {
    for (i=0; i < no_oop_types; i++) {
      for (j=0; j < 3; j++) ooptypes[i].params[j] = 0.0;
      for (j=0; j < 4; j++)
        strncpy(potential_types[j],
                atomtypes[ooptypes[i].types[j]].potential,5);

      k = find_improper_body_data(potential_types,ff_oop,&rearrange);
      if (k < 0) {
        get_equivs(5,potential_types,equiv_types);

        if (pflag > 2) {
          printf(" Using equivalences for oop %s %s %s %s -> %s %s %s %s\n",
                 potential_types[0],potential_types[1],
                 potential_types[2],potential_types[3],
                 equiv_types[0],equiv_types[1],
                 equiv_types[2],equiv_types[3]);
        }
        k = find_improper_body_data(equiv_types,ff_oop,&rearrange);
      }
      if (k < 0) {
        printf(" Unable to find oop data for %s %s %s %s\n",
               potential_types[0],
               potential_types[1],potential_types[2],potential_types[3]);
        condexit(22);
      } else {
        ooptypes[i].params[0] = ff_oop.data[k].ff_param[0];
        if (ff_oop.data[k].ff_param[2] == 0.0)
          ooptypes[i].params[1] = 1.0;
        else if (ff_oop.data[k].ff_param[2] == 180.0)
          ooptypes[i].params[1] = -1.0;
        else {
          printf(" Non planar phi0 for %s %s %s %s\n",
                 potential_types[0],potential_types[1],
                 potential_types[2],potential_types[3]);
          ooptypes[i].params[1] = 0.0;
        }
        ooptypes[i].params[2] = ff_oop.data[k].ff_param[1];
        if (rearrange > 0) rearrange_improper(i,rearrange);
      }
    }
  }

  if (forcefield & FF_TYPE_CLASS2) {
    for (i=0; i < no_oop_types; i++) {
      for (j=0; j < 3; j++)
        ooptypes[i].params[j] = 0.0;
      for (j=0; j < 4; j++)
        strncpy(potential_types[j],
                atomtypes[ooptypes[i].types[j]].potential,5);
      k = find_trigonal_body_data(potential_types,ff_oop);
      if (k < 0) {
        get_equivs(5,potential_types,equiv_types);
        if (pflag > 2) {
          printf(" Using equivalences for oop %s %s %s %s -> %s %s %s %s\n",
                 potential_types[0],potential_types[1],
                 potential_types[2],potential_types[3],
                 equiv_types[0],equiv_types[1],
                 equiv_types[2],equiv_types[3]);
        }
        k = find_trigonal_body_data(equiv_types,ff_oop);
      }
      if (k < 0) {
        printf(" Unable to find oop data for %s %s %s %s\n",
               potential_types[0],
               potential_types[1],potential_types[2],potential_types[3]);
        condexit(23);
      } else {
        for (j=0; j < 2; j++)
          ooptypes[i].params[j] = ff_oop.data[k].ff_param[j];
      }
    }
  }

  if (pflag > 2) {
    printf("\n OOP Types and  Parameters\n");
    for (i=0; i < no_oop_types; i++) {
      for (j=0; j < 4; j++)
        printf(" %-3s",atomtypes[ooptypes[i].types[j]].potential);
      for (j=0; j < 3; j++)
        printf(" %8.4f",ooptypes[i].params[j]);
      printf("\n");
    }
  }


  /**********************************************************************/
  /*                                                                    */
  /*   Find parameters for angleangle types (Class II only)             */
  /*                                                                    */
  /*   This is somewhat complicated in that one set of four types       */
  /*   a b c d has three angleangle combinations so for each type       */
  /*   the program needs to find three sets of parameters by            */
  /*   progressively looking for data for different permutations of     */
  /*   a c and d                                                        */
  /*                                                                    */
  /**********************************************************************/

  if (forcefield & FF_TYPE_CLASS2) {

    for (i=0; i < no_oop_types; i++) {

      for (j=0; j < 6; j++) ooptypes[i].angleangle_params[j] = 0.0;

      for (j=0; j < 4; j++)
        strncpy(potential_types[j],
                atomtypes[ooptypes[i].types[j]].potential,5);


      tabc = get_t0(ooptypes[i].types[0],
                    ooptypes[i].types[1],
                    ooptypes[i].types[2]);

      tabd = get_t0(ooptypes[i].types[0],
                    ooptypes[i].types[1],
                    ooptypes[i].types[3]);
      tcbd = get_t0(ooptypes[i].types[2],
                    ooptypes[i].types[1],

                    ooptypes[i].types[3]);

      ooptypes[i].angleangle_params[3] = tabc;
      ooptypes[i].angleangle_params[4] = tcbd;
      ooptypes[i].angleangle_params[5] = tabd;

      k = find_angleangle_data(potential_types,ff_angang,kloc);
      if (k < 0) {
        get_equivs(5,potential_types,equiv_types);
        if (pflag > 2) {
          printf(" Using equivalences for angleangle %s %s %s %s -> %s %s %s %s\n",
                 potential_types[0],potential_types[1],
                 potential_types[2],potential_types[3],
                 equiv_types[0],equiv_types[1],
                 equiv_types[2],equiv_types[3]);
          k = find_angleangle_data(equiv_types,ff_angang,kloc);
        }
      }
      if (k < 0) {
        printf(" Unable to find angleangle data for %s %s %s %s\n",
               potential_types[0],
               potential_types[1],potential_types[2],potential_types[3]);
        condexit(24);
      } else {
        for (j=0; j < 3; j++) {
          if (kloc[j] > -1)
            ooptypes[i].angleangle_params[j] = ff_angang.data[kloc[j]].ff_param[0];
        }
      }
    }

    for (i=0; i < no_angleangle_types; i++) {
      for (j=0; j < 6; j++) angleangletypes[i].params[j] = 0.0;
      for (j=0; j < 4; j++)
        strncpy(potential_types[j],
                atomtypes[angleangletypes[i].types[j]].potential,5);

      tabc = get_t0(angleangletypes[i].types[0],
                    angleangletypes[i].types[1],
                    angleangletypes[i].types[2]);
      tabd = get_t0(angleangletypes[i].types[0],
                    angleangletypes[i].types[1],
                    angleangletypes[i].types[3]);
      tcbd = get_t0(angleangletypes[i].types[2],
                    angleangletypes[i].types[1],
                    angleangletypes[i].types[3]);

      angleangletypes[i].params[3] = tabc;
      angleangletypes[i].params[4] = tcbd;
      angleangletypes[i].params[5] = tabd;

      k = find_angleangle_data(potential_types,ff_angang,kloc);
      if (k < 0) {
        get_equivs(5,potential_types,equiv_types);
        if (pflag > 2) {
          printf("Using equivalences for angleangle %s %s %s %s -> %s %s %s %s\n",
                 potential_types[0],potential_types[1],
                 potential_types[2],potential_types[3],
                 equiv_types[0],equiv_types[1],
                 equiv_types[2],equiv_types[3]);
        }
        k = find_angleangle_data(equiv_types,ff_angang,kloc);
      }
      if (k < 0) {
        printf(" Unable to find angleangle data for %s %s %s %s\n",
               potential_types[0],
               potential_types[1],potential_types[2],potential_types[3]);
        condexit(25);
      } else {
        for (j=0; j < 3; j++) {
          if (kloc[j] > -1)
            angleangletypes[i].params[j] =
              ff_angang.data[kloc[j]].ff_param[0];
        }
      }
    }
    if (pflag > 2) {
      printf("\n AngleAngle Types and  Parameters\n");
      for (i=0; i < no_oop_types; i++) {
        for (j=0; j < 4; j++)
          printf(" %-3s",atomtypes[ooptypes[i].types[j]].potential);
        for (j=0; j < 6; j++)
          printf(" %8.4f",ooptypes[i].angleangle_params[j]);
        printf("\n");
      }
      for (i=0; i < no_angleangle_types; i++) {
        for (j=0; j < 4; j++)
          printf(" %-3s",atomtypes[angleangletypes[i].types[j]].potential);
        for (j=0; j < 6; j++) printf(" %8.4f",angleangletypes[i].params[j]);
        printf("\n");
      }
    }
  }
}

int find_improper_body_data(char types1[][5],struct FrcFieldItem item,
                            int *rearrange_ptr)
{
  int k,backwards;
  char mirror_types[4][5];

  backwards = 0;

  /* a b c d */

  *rearrange_ptr = 0;
  k = find_match(4,types1,item,&backwards);
  if (k >= 0) return k;

  /* a b d c */

  *rearrange_ptr = 1;
  strncpy(mirror_types[0],types1[0],5);
  strncpy(mirror_types[1],types1[1],5);
  strncpy(mirror_types[2],types1[3],5);
  strncpy(mirror_types[3],types1[2],5);
  k = find_match(4,mirror_types,item,&backwards);
  if (k >= 0) return k;

  /* d b a c */

  *rearrange_ptr = 2;
  strncpy(mirror_types[0],types1[3],5);
  strncpy(mirror_types[2],types1[0],5);
  strncpy(mirror_types[3],types1[2],5);
  k = find_match(4,mirror_types,item,&backwards);
  if (k >= 0) return k;

  /* d b c a */

  *rearrange_ptr = 3;
  strncpy(mirror_types[2],types1[2],5);
  strncpy(mirror_types[3],types1[0],5);
  k = find_match(4,mirror_types,item,&backwards);
  if (k >= 0) return k;

  /* c b a d */

  *rearrange_ptr = 4;
  strncpy(mirror_types[0],types1[2],5);
  strncpy(mirror_types[2],types1[0],5);
  strncpy(mirror_types[3],types1[3],5);
  k = find_match(4,mirror_types,item,&backwards);
  if (k >= 0) return k;

  /* c b d a */

  *rearrange_ptr = 5;
  strncpy(mirror_types[2],types1[3],5);
  strncpy(mirror_types[3],types1[0],5);
  k = find_match(4,mirror_types,item,&backwards);
  return k;
}

void rearrange_improper(int ooptype,int rearrange)
{
  int i,j,temp[4];

  for (i=0; i < 4; i++) temp[i] = ooptypes[ooptype].types[i];

  switch (rearrange) {
  case 1:
    ooptypes[ooptype].types[0] = temp[0];
    ooptypes[ooptype].types[2] = temp[3];
    ooptypes[ooptype].types[3] = temp[2];
    for (i=0; i < total_no_oops; i++) {
      if (oops[i].type == ooptype) {
        for (j=0; j < 4; j++) temp[j] = oops[i].members[j];
        oops[i].members[2] = temp[3];
        oops[i].members[3] = temp[2];
      }
    }
    break;
  case 2:
    ooptypes[ooptype].types[0] = temp[3];
    ooptypes[ooptype].types[2] = temp[0];
    ooptypes[ooptype].types[3] = temp[2];
    for (i=0; i < total_no_oops; i++) {
      if (oops[i].type == ooptype) {
        for (j=0; j < 4; j++) temp[j] = oops[i].members[j];
        oops[i].members[0] = temp[3];
        oops[i].members[2] = temp[0];
        oops[i].members[3] = temp[2];
      }
    }
    break;
  case 3:
    ooptypes[ooptype].types[0] = temp[3];
    ooptypes[ooptype].types[2] = temp[2];
    ooptypes[ooptype].types[3] = temp[0];
    for (i=0; i < total_no_oops; i++) {
      if (oops[i].type == ooptype) {
        for (j=0; j < 4; j++) temp[j] = oops[i].members[j];
        oops[i].members[0] = temp[3];
        oops[i].members[2] = temp[2];
        oops[i].members[3] = temp[0];
      }
    }
    break;
  case 4:
    ooptypes[ooptype].types[0] = temp[2];
    ooptypes[ooptype].types[2] = temp[0];
    ooptypes[ooptype].types[3] = temp[3];
    for (i=0; i < total_no_oops; i++) {
      if (oops[i].type == ooptype) {
        for (j=0; j < 4; j++) temp[j] = oops[i].members[j];
        oops[i].members[0] = temp[2];
        oops[i].members[2] = temp[0];
        oops[i].members[3] = temp[3];
      }
    }
    break;
  case 5:
    ooptypes[ooptype].types[0] = temp[2];
    ooptypes[ooptype].types[2] = temp[3];
    ooptypes[ooptype].types[3] = temp[0];
    for (i=0; i < total_no_oops; i++) {
      if (oops[i].type == ooptype) {
        for (j=0; j < 4; j++) temp[j] = oops[i].members[j];
        oops[i].members[0] = temp[2];
        oops[i].members[2] = temp[3];
        oops[i].members[3] = temp[0];
      }
    }
    break;
  default:
    break;
  }
}

int find_trigonal_body_data(char types1[][5],struct FrcFieldItem item)
{
  int k,backwards;
  char mirror_types[4][5];

  backwards = -1;

  /* a b c d */

  k = find_match(4,types1,item,&backwards);
  if (k >= 0) return k;

  /* a b d c */

  strncpy(mirror_types[0],types1[0],5);
  strncpy(mirror_types[1],types1[1],5);
  strncpy(mirror_types[2],types1[3],5);
  strncpy(mirror_types[3],types1[2],5);
  k = find_match(4,mirror_types,item,&backwards);
  if (k >= 0) return k;

  /* d b a c */

  strncpy(mirror_types[0],types1[3],5);
  strncpy(mirror_types[2],types1[0],5);
  strncpy(mirror_types[3],types1[2],5);
  k = find_match(4,mirror_types,item,&backwards);
  if (k >= 0) return k;

  /* d b c a */

  strncpy(mirror_types[2],types1[2],5);
  strncpy(mirror_types[3],types1[0],5);
  k = find_match(4,mirror_types,item,&backwards);
  if (k >= 0) return k;
  /* c b a d */

  strncpy(mirror_types[0],types1[2],5);
  strncpy(mirror_types[2],types1[0],5);
  strncpy(mirror_types[3],types1[3],5);
  k = find_match(4,mirror_types,item,&backwards);
  if (k >= 0) return k;

  /* c b d a */

  strncpy(mirror_types[2],types1[3],5);
  strncpy(mirror_types[3],types1[0],5);
  k = find_match(4,mirror_types,item,&backwards);
  return k;
}

int find_angleangle_data(char types1[][5],struct FrcFieldItem item,int kloc[3])
{
  int k,backwards = -1;
  char mirror_types[4][5];

  strncpy(mirror_types[1],types1[1],5);

  /* go for first parameter a b c d or d b c a */

  k = find_match(4,types1,item,&backwards);
  if (k < 0) {
    strncpy(mirror_types[0],types1[3],5);
    strncpy(mirror_types[2],types1[2],5);
    strncpy(mirror_types[3],types1[0],5);
    k = find_match(4,mirror_types,item,&backwards);
  }
  kloc[0] = k;

  /* go for second parameter d b a c or c b a d */

  strncpy(mirror_types[0],types1[3],5);
  strncpy(mirror_types[2],types1[0],5);
  strncpy(mirror_types[3],types1[2],5);
  k = find_match(4,mirror_types,item,&backwards);
  if (k < 0) {
    strncpy(mirror_types[0],types1[2],5);
    strncpy(mirror_types[3],types1[3],5);
    k = find_match(4,mirror_types,item,&backwards);
  }
  kloc[1] = k;

  /* go for third parameter a b d c or c b d a */

  strncpy(mirror_types[0],types1[0],5);
  strncpy(mirror_types[2],types1[3],5);
  strncpy(mirror_types[3],types1[2],5);
  k = find_match(4,mirror_types,item,&backwards);
  if (k < 0) {
    strncpy(mirror_types[0],types1[2],5);
    strncpy(mirror_types[3],types1[0],5);
    k = find_match(4,mirror_types,item,&backwards);
  }
  kloc[2] = k;
  k = 0;
  if ((kloc[0] < 0) && (kloc[1] < 0) && (kloc[2] < 0)) k = -1;
  return k;
}

int find_match(int n, char types1[][5],struct FrcFieldItem item,int
               *backwards_ptr)
{
  int k,match;

  match = 0;
  k=0;

  /* Try for an exact match (no wildcards) first */

  while (!match && (k < item.entries)) {
    if (match_types(n, 0,types1,item.data[k].ff_types,backwards_ptr) == 1)
      match = 1;
    else
      k++;
  }

  /* Try again - allow wildcard matching  */

  if (!match) {
    k=0;
    while (!match && (k < item.entries)) {
      if (match_types(n,1,types1,item.data[k].ff_types,backwards_ptr) == 1)
        match = 1;
      else
        k++;
    }
  }
  if (match) return k;
  else return -1;
}

int match_types(int n,int wildcard,char types1[][5],char types2[][5],
                int *backwards_ptr)
{
  int k,match;

  /* Routine to match short arrays of characters strings which contain
     atom potential types. The arrays range from 1 to 4 (VDW or equivalences,
     bond, angle, dihedrals or oops). There are potentially four ways the
     arrays can match: exact match (forwards), exact match when one array is
     run backwards (backwards), forwards with wildcard character match allowed
     (forwards *) and finally backwards with wildcard character match
     (backwards *). If the variable, backwards (pointed by backwards_ptr)
     is -1, then the backwards options are not to be used (such when
     matching oop types)
  */


  if (wildcard == 0) {

  /* forwards */

    k=0;
    match = 1;
    while (match && (k < n)) {
      if (strncmp(types1[k],types2[k],5) == 0)
        k++;
      else
        match = 0;
    }
  } else {

  /* forwards * */

    k=0;

    match = 1;
    while (match && (k < n)) {
      if ((strncmp(types1[k],types2[k],5) == 0) ||
          (types2[k][0] == '*'))
        k++;
      else
        match = 0;
    }
  }

  if (match) {
    *backwards_ptr = 0;
    return 1;
  }
  if ((n < 2) || (*backwards_ptr == -1)) return 0;

  if (wildcard == 0) {

  /* backwards */

    k=0;
    match = 1;
    while (match && (k < n)) {
      if (strncmp(types1[n-k-1],types2[k],5) == 0)
        k++;
      else
        match = 0;
    }
  } else {

  /* backwards * */

    k=0;
    match = 1;
    while (match && (k < n)) {
      if ((strncmp(types1[n-k-1],types2[k],5) == 0) ||
          (types2[k][0] == '*')                   )
        k++;
      else
        match = 0;
    }
  }

  if (match) {
    *backwards_ptr = 1;
    return 1;
  } else return 0;
}

double get_r0(int typei,int typej)
{
  int k,match;
  double r;

  k=0;
  match=0;
  r = 0.0;

  while (!match && (k < no_bond_types)) {
    if (((typei == bondtypes[k].types[0]) &&
         (typej == bondtypes[k].types[1])) ||
        ((typej == bondtypes[k].types[0]) &&
         (typei == bondtypes[k].types[1]))   ) {
      r = bondtypes[k].params[0];
      match = 1;
    } else k++;
  }

  if (match == 0)
    printf(" Unable to find r0 for types %d %d\n",typei,typej);
  return r;
}

double get_t0(int typei,int typej,int typek)
{
  int k,match;
  double theta;

  k=0;
  match=0;
  theta = 0.0;

  while (!match && (k < no_angle_types)) {
    if (((typei == angletypes[k].types[0]) &&
         (typej == angletypes[k].types[1]) &&
         (typek == angletypes[k].types[2])) ||
        ((typek == angletypes[k].types[0]) &&
         (typej == angletypes[k].types[1]) &&
         (typei == angletypes[k].types[2]))   ) {
      theta = angletypes[k].params[0];
      match = 1;
    } else k++;
  }

  if (match == 0)
    printf(" Unable to find t0 for types %d %d %d\n",
           typei,typej,typek);
  return theta;
}

int quo_cp()
{
  char cp[] = "cp  ";
  int i,type,found;

  i = 0;
  type = -1;
  found = 0;

  while (!found && (i < no_atom_types)) {
    if (strncmp(atomtypes[i].potential,cp,2) == 0) {
      found = 1;
      type = i;
    } else i++;
  }

  return type;
}

void get_equivs(int ic,char potential_types[][5],char equiv_types[][5])
{
  int i,k;
  switch (ic) {
  case 1:
    k = find_equiv_type(potential_types[0]);
    if (k > -1) strncpy(equiv_types[0],equivalence.data[k].ff_types[1],5);
    break;

  case 2:
    for (i=0; i < 2; i++) {
      k = find_equiv_type(potential_types[i]);
      if (k > -1) strncpy(equiv_types[i],equivalence.data[k].ff_types[2],5);
    }
    break;
  case 3:
    for (i=0; i < 3; i++) {
      k = find_equiv_type(potential_types[i]);
      if (k > -1) strncpy(equiv_types[i],equivalence.data[k].ff_types[3],5);
    }
    break;
  case 4:
    for (i=0; i < 4; i++) {
      k = find_equiv_type(potential_types[i]);
      if (k > -1) strncpy(equiv_types[i],equivalence.data[k].ff_types[4],5);
    }
    break;

  case 5:
    for (i=0; i < 4; i++) {
      k = find_equiv_type(potential_types[i]);
      if (k > -1)
        strncpy(equiv_types[i],equivalence.data[k].ff_types[5],5);
    }
    break;
  default:
    printf(" Requesting equivalences of unsupported type: %d\n",ic);
    condexit(26);
    break;
  }
  return;
}

int find_equiv_type(char potential_type[5])
{
  int j,k,match;

  j = -1;
  k = 0;
  match = 0;

  while (!match && (k < equivalence.entries)) {
    if (strncmp(potential_type,
                equivalence.data[k].ff_types[0],5) == 0) {
      match = 1;
      j = k;
    } else {
      k++;
    }
  }
  if (j < 0)
    printf(" Unable to find equivalent type for %s\n",potential_type);
  return j;
}
