label_flex.c

来自「最经典的分子对结软件」· C语言 代码 · 共 1,658 行 · 第 1/3 页

          "INFLEXIBLE"
      );
    }

    if (label_flex->drive_flag)
      fprintf (global.outfile,
        "\nTotal conformations (ignoring symmetry and clashes) for %s %.1g\n",
        molecule->info.name, conf_total);
  }
}
 

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

Subroutine to identify all bonds in rings.

11/96 te

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

void detect_rings (MOLECULE *molecule)
{
  int i, j;

  int next_atom_level (int, int, int, MOLECULE *);

/*
* Reset atom flags and bond ring flags
* 4/97 te
*/
  for (i = 0; i < molecule->total.atoms; i++)
    molecule->atom[i].flag = 0;

  for (i = 0; i < molecule->total.bonds; i++)
    molecule->bond[i].ring_flag = 0;

/*
* Identify ring bonds
* 4/97 te
*/
  next_atom_level (1, 0, 0, molecule);

/*
* Also identify bonds specified in RIGID set, if present
* 4/97 te
*/
  for (i = 0; i < molecule->total.sets; i++)
  {
    if ((molecule->set[i].name != NULL) &&
      !strcmp (molecule->set[i].name, "RIGID"))
    {
      if ((molecule->set[i].type != NULL) &&
        !strcmp (molecule->set[i].type, "STATIC"))
      {
        if ((molecule->set[i].obj_type != NULL) &&
          !strcmp (molecule->set[i].obj_type, "BONDS"))
        {
          if (molecule->set[i].member_total > 0)
          {
            for (j = 0; j < molecule->set[i].member_total; j++)
            {
              if ((molecule->set[i].member[j] >= 0) &&
                (molecule->set[i].member[j] < molecule->total.bonds))
                molecule->bond[molecule->set[i].member[j]].ring_flag = TRUE;

              else
                exit (fprintf (global.outfile,
                  "ERROR detect_rings: RIGID set contains invalid bonds\n"));
            }
          }

          else
            exit (fprintf (global.outfile,
              "ERROR detect_rings: RIGID set empty\n"));
        }

        else
          exit (fprintf (global.outfile,
            "ERROR detect_rings: RIGID set must be BONDS obj_type\n"));
      }

      else
        exit (fprintf (global.outfile,
          "ERROR detect_rings: RIGID set must be STATIC type\n"));
    }
  }

/*
  for (i = 0; i < molecule->total.bonds; i++)
    if (molecule->bond[i].ring_flag)
      fprintf (global.outfile, "Ring bond %d %s %s\n", i + 1,
        molecule->atom[molecule->bond[i].origin].name,
        molecule->atom[molecule->bond[i].target].name);
*/
}

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

Recursive subroutine that traverses all atom linkage paths in a molecule.
The level of recursion is recorded for each atom as it is traversed.
The recursion level is reported back, unless an atom is re-encountered
during a higher level of recursion; then, the level of the previously seen
atom is reported.

A ring bond is identified when an atom is linked to an atom previously seen.

4/97 te

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

int next_atom_level
(
  int level,
  int current_atom,
  int previous_atom,
  MOLECULE *molecule
)
{
  int i;
  int bond_id;			/* Bond id linking current neighbor */
  int neighbor_flag;		/* Flag of neighbor atom */

/*
* If the current atom hasn't been flagged, then check the flags
* of its neighbors
* 4/97 te
*/
  if (molecule->atom[current_atom].flag == 0)
  {
    molecule->atom[current_atom].flag = level;        

    for (i = 0; i < molecule->atom[current_atom].neighbor_total; i++)
    {
      if (molecule->atom[current_atom].neighbor[i].id != previous_atom) 
      {
        neighbor_flag =
          next_atom_level
          (
            level + 1,
            molecule->atom[current_atom].neighbor[i].id,
            current_atom,
            molecule
          );

        if (neighbor_flag <= level)
        {
          bond_id = molecule->atom[current_atom].neighbor[i].bond_id;
          molecule->bond[bond_id].ring_flag = TRUE;
        }

        if (neighbor_flag < molecule->atom[current_atom].flag)
          molecule->atom[current_atom].flag = neighbor_flag;
      }
    }
  }

  return molecule->atom[current_atom].flag;
}


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

int check_peripheral_torsion
(
  MOLECULE	*molecule,
  int		torsion_id
)
{
  int atom_id;

  atom_id = molecule->torsion[torsion_id].origin_neighbor;
  if (molecule->atom[atom_id].heavy_flag == FALSE)
    return TRUE;

  atom_id = molecule->torsion[torsion_id].target_neighbor;
  if (molecule->atom[atom_id].heavy_flag == FALSE)
    return TRUE;

  return FALSE;
}


/* /////////////////////////////////////////////////////////////////// */
/*
Routine to assign torsions to segments and segments to layers.
It also loops through all anchor fragments.

Return values:
  TRUE	successful identification of anchor segment
  EOF	unable to identify any more anchor segments

12/96 te

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

int get_anchor
(
  LABEL_FLEX	*label_flex,
  MOLECULE	*mol_init,
  MOLECULE	*mol_anch,
  int		anchor
) 
{
/*
* Initialize variables
* 12/96 te
*/
  if (anchor == 0)
  {
/*
*   Assign flexible labels
*   12/96 te
*/
    if (label_flex->flag)
    {
      mol_init->transform.tors_flag = FALSE;
      assign_flex_labels (label_flex, mol_init);

      if (mol_init->total.torsions > label_flex->max_torsions)
      {
        if (global.output_volume != 't')
          fprintf(global.outfile, "Skipped %s (%d rotatable bonds).\n",
            mol_init->info.name, mol_init->total.torsions);

        return EOF;
      }
    }

    get_segments (label_flex, mol_init);
    copy_molecule (mol_anch, mol_init);
  }

/*
* Exit if return visit, but multiple anchors not requested
* 12/96 te
*/
  else if (!label_flex->multiple_anchors)
    return EOF;

  else
    copy_segments (mol_anch, mol_init);

/*
* Find next anchor
* 12/96 te
*/
  if (get_anchor_segment (label_flex, mol_anch, anchor) != TRUE)
    return EOF;

  initialize_segments (label_flex, mol_anch);
  initialize_anchor_layer (label_flex, mol_anch);

  return TRUE;
}


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

Routine to divide a molecule into rigid segments.

11/96 te

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

void get_segments (LABEL_FLEX *label_flex, MOLECULE *mol_init)
{
  int i, j;
  int atom_id;			/* Atom id */
  int bond_id;			/* Bond id */
  int torsion_id;		/* Torsion id */
  int origin;			/* Origin atom id */
  int target;			/* Target atom id */
  int segment;			/* Current segment id */
  int neighbor;			/* Neighboring segment id */
  int neighbor_id;		/* Current position in neighbor list */

  static SEARCH search = {0};

/*
* Allocate more than enough space for segments
* 11/96 te
*/
  mol_init->total.segments = mol_init->total.torsions + 1;
  reallocate_segments (mol_init);

  for (i = 0; i < mol_init->total.segments; i++)
  {
    mol_init->segment[i].atom_total = mol_init->total.atoms;
    mol_init->segment[i].neighbor_total = mol_init->total.torsions;

    reallocate_segment_atoms (&mol_init->segment[i]);
    reallocate_segment_neighbors (&mol_init->segment[i]);

    reset_segment (&mol_init->segment[i]);
  }

/*
* Identify rigid segments separated by rotatable bonds
* 11/96 te
*/
  if (mol_init->total.torsions == 0)
  {
    for (i = 0; i < mol_init->total.atoms; i++)
    {
      mol_init->segment[0].atom[i] = i;
      mol_init->atom[i].segment_id = 0;
    }

    mol_init->segment[0].atom_total = mol_init->total.atoms;
  }

  else
  {
    for
    (
      i = mol_init->total.segments = 0;
      breadth_search
      (
        &search,
        mol_init->atom,
        mol_init->total.atoms,
        get_atom_neighbor,
        NULL,
        &i, 1, NEITHER, i
      ) != EOF;
      i++
    )
    {
      origin = get_search_origin (&search, NEITHER, NEITHER, NEITHER);
      target = get_search_target (&search, NEITHER, NEITHER, NEITHER);

/*
*     If this is the first atom, then initialize segment list only
*     11/96 te
*/
      if (i == 0)
      {
        mol_init->atom[target].segment_id = 0;
        mol_init->segment[0].atom[0] = target;
        mol_init->segment[0].atom_total = 1;
        mol_init->total.segments = 1;
        continue;
      }

/*
*     Check if this linkage is a rotatable bond
*     11/96 te
*/
      for
      (
        j = 0, torsion_id = NEITHER;
        (j < mol_init->total.torsions) && (torsion_id == NEITHER);
        j++
      )
        if (mol_init->torsion[j].flex_id > 0)
        {
          bond_id = mol_init->torsion[j].bond_id;

          if (((mol_init->bond[bond_id].origin == origin) &&
            (mol_init->bond[bond_id].target == target)) ||
            ((mol_init->bond[bond_id].origin == target) &&
            (mol_init->bond[bond_id].target == origin)))
            torsion_id = j;
        }

/*
*     If the link is flexible, then increment the segment total,
*     update the segment neighbor lists and atom list, and atom segment id
*     11/96 te
*/
      if (torsion_id != NEITHER)
      {
        segment = mol_init->total.segments++;
        neighbor = mol_init->atom[origin].segment_id;

        neighbor_id = mol_init->segment[segment].neighbor_total;
        mol_init->segment[segment].neighbor[neighbor_id].id = neighbor;
        mol_init->segment[segment].neighbor_total++;

        neighbor_id = mol_init->segment[neighbor].neighbor_total;
        mol_init->segment[neighbor].neighbor[neighbor_id].id = segment;
        mol_init->segment[neighbor].neighbor_total++;

        atom_id = mol_init->segment[segment].atom_total;
        mol_init->segment[segment].atom[atom_id] = target;
        mol_init->segment[segment].atom_total++;

        mol_init->atom[target].segment_id = segment;

        mol_init->segment[segment].periph_flag =
          mol_init->torsion[torsion_id].periph_flag;
      }

/*
*     Otherwise, update the segment atom list and atom segment id
*     11/96 te
*/
      else
      {
        segment = mol_init->atom[origin].segment_id;

        atom_id = mol_init->segment[segment].atom_total;
        mol_init->segment[segment].atom[atom_id] = target;
        mol_init->segment[segment].atom_total++;

        mol_init->atom[target].segment_id = segment;
      }
    }
  }

/*
* Determine the size of each segment
* 1/97 te
*/
  for (i = 0; i < mol_init->total.segments; i++)
  {
    mol_init->segment[i].heavy_total = 0;

    for (j = 0; j < mol_init->segment[i].atom_total; j++)
    {
      atom_id = mol_init->segment[i].atom[j];

      if (mol_init->atom[atom_id].heavy_flag == TRUE)
        mol_init->segment[i].heavy_total++;
    }
  }

  if (label_flex->flag && (global.output_volume == 'v'))
  {
    fprintf (global.outfile, "\nInitial segments\n");

    for (i = 0; i < mol_init->total.segments; i++)
    {
      fprintf (global.outfile, "  segment %d\n", i + 1);
      fprintf (global.outfile, "    neighbors:");

      for (j = 0; j < mol_init->segment[i].neighbor_total; j++)
        fprintf (global.outfile, " %d",
          mol_init->segment[i].neighbor[j].id + 1);

      fprintf (global.outfile, "\n    atoms    :");

      for (j = 0; j < mol_init->segment[i].atom_total; j++)
        fprintf (global.outfile, " %s",
          mol_init->atom[mol_init->segment[i].atom[j]].name);

      fprintf (global.outfile, "\n");
    }
  }
}


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

Routine to identify segments to use as anchors.

Return values:
	TRUE:	next anchor found
	EOF:	unable to find any more anchors

12/96 te

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

int get_anchor_segment
(
  LABEL_FLEX	*label_flex,
  MOLECULE	*molecule,
  int		anchor_count
)
{
  int i;
  int size;
  int anchor_found = FALSE;	/* Flag for whether anchor found */

  if (anchor_count > 0)
  {
    if (label_flex->multiple_anchors == FALSE)
      return EOF;
  }

  else
  {
    if (label_flex->anchor_flag == FALSE)
      return TRUE;

    else if (molecule->transform.anchor_atom != NEITHER)
      exit (fprintf (global.outfile,
        "ERROR get_anchor_segment: Molecule %s %s\n"
        "  This molecule has an anchor atom and transformation from input.\n"
        "  Convert input file to non-PTR format before processing.\n",
        molecule->info.name, molecule->info.comment));
  }

/*
* Check if an anchor was specified in input
* 1/97 te
*/
  for (i = 0; i < molecule->total.sets; i++)
    if (!strcmp (molecule->set[i].name, "ANCHOR"))
    {
      if ((strcmp (molecule->set[i].type, "STATIC")) ||
        (strcmp (molecule->set[i].obj_type, "ATOMS")))
        exit (fprintf (global.outfile,
          "ERROR get_anchor_segment: Molecule %s %s\n"
          "  This molecule has an ANCHOR set in input file.\n"
          "  ANCHOR set must be of type STATIC ATOMS.\n",
          molecule->info.name, molecule->info.comment));

      if (label_flex->multiple_anchors == TRUE)
        exit (fprintf (global.outfile,
          "ERROR get_anchor_segment: Molecule %s %s\n"
          "  This molecule has an ANCHOR atom set in input file.\n"
          "  Delete ANCHOR set before processing with multiple_anchors.\n",
          molecule->info.name, molecule->info.comment));

      if
      (
        (molecule->set[i].member_total > 0) &&
        (molecule->set[i].member[0] >= 0) &&
        (molecule->set[i].member[0] < molecule->total.atoms)
      )
      {
        molecule->transform.anchor_atom = molecule->set[i].member[0];
        return TRUE;
      }

      else
        exit (fprintf (global.outfile,
          "ERROR get_anchor_segment: Molecule %s %s\n"
          "  This molecule has an ANCHOR atom set in input file.\n"
          "  The ANCHOR set contains no atoms or incorrect atoms.\n",
          molecule->info.name, molecule->info.comment));
    }

/*
* If multiple anchors allowed, then find the next segment
* that satisfies the anchor size cutoff
* 1/97 te
*/
  if ((anchor_found == FALSE) &&
    (label_flex->multiple_anchors == TRUE))
  {
    for
    (
      anchor_found = FALSE,
        i = (anchor_count > 0) ?
          molecule->atom[molecule->transform.anchor_atom].segment_id + 1 : 0;
      (anchor_found == FALSE) &&
        (i < molecule->total.segments);
      i++