match.c

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

int match_distance (MATCH *match)
{
  int i, j;			/* Counter variable */
  EDGE edge;			/* Edge under consideration */
  EDGE iedge, jedge;		/* Edge comparisons */
  int valid_clique;		/* Flag for whether an valid_clique found */

  int match_chirality (MATCH *, int);

/*
* If a prior clique existed, then erase the last node and
* decrement the clique size
* 11/96 te
*/
  if (match->clique.edge_total > 0)
  {
    match->clique.edge_total--;
    match->clique.edge[match->clique.edge_total].node = 0;
    match->clique.edge[match->clique.edge_total].residual = 0;
  }

/*
* Perform clique detection until valid clique is found
* 11/96 te
*/
  for (valid_clique = FALSE; !valid_clique;)
  {
/*
*   Can we expand the current clique?  Check if:
*   1. The clique is smaller than the maximum size
*   2. There are unused nodes in the clique adjacency list
*   11/96 te
*/
    if
    (
      (match->clique.edge_total < match->clique_size_max) &&
      (match->adjacent[match->clique.edge_total] <
        match->clique_adjacency_total[match->clique.edge_total])
    )
    {
/*
*     Identify the next adjacent node and increment the adjacency list
*     11/96 te
*/
      edge = match->clique_adjacency_list
        [match->clique.edge_total]
        [match->adjacent[match->clique.edge_total]++];

/*
*     STAGE 1: The node is geometrically feasible, if
*     1. The node passes the critical cluster filter:
*        a. The filter is not requested, OR
*        b. The filter indicates valid node, AND
*     2. The clique has proper chirality:
*        a. Chirality checking is not necessary, OR
*        b. The node is not being added to a three-node clique, OR
*        c. The chirality of the new clique is correct, OR
*        d. The ligand is allowed to be reflected
*     11/96 te
*/
      if
      (
        (
          !match->critical_flag ||
          (((match->clique.edge_total == 0) &&
            match->critical_filter[0][edge.node]) ||
          ((match->clique.edge_total > 0) &&
            match->critical_filter
              [match->receptor_site.atom
                [match->clique.edge[match->clique.edge_total - 1].node /
                  match->ligand_center.total.atoms].subst_id + 1][edge.node]))
        ) &&
        (
          !match->chiral_flag ||
          (match->clique.edge_total != 3) ||
          !(match->clique.reflect_flag = match_chirality (match, edge.node)) ||
          match->reflect_flag
        )
      )
      {
/*
*       STAGE 2: The node is worth adding, if:
*       1. The current node isn't forming a redundant clique:
*          a. We don't care, OR
*          b. It the first in the clique, OR
*          c. The node is numerically above the previously added node, AND
*       2. The current node passes the chemical filter
*          a. The filter is not requested, OR
*          a. The filter indicates a valid node
*       11/96 te
*/
        if
        (
          (
            !match->degeneracy_flag ||
            !match->clique.edge_total ||
            (edge.node > match->clique.edge[match->clique.edge_total - 1].node)
          ) &&
          (
            !match->chemical_flag ||
            match->chemical_filter[edge.node]
          )
        )
        {
/*
*         Record the node in the clique and increment the clique size.
*         Also, document the finding of a clique of this size,
*         and erase the record of the next larger one.
*         11/96 te
*/
          match->clique.edge[match->clique.edge_total++] = edge;

          if (match->degeneracy_flag)
          {
            match->degenerate[match->clique.edge_total - 1] = TRUE;
            match->degenerate[match->clique.edge_total] = FALSE;
          }

/*
*         Construct an adjacency list for the new clique
*         11/96 te
*/
          for
          (
            i = j = match->adjacent[match->clique.edge_total] = 
              match->clique_adjacency_total[match->clique.edge_total] = 0;
            (i < match->clique_adjacency_total[match->clique.edge_total - 1]) &&
              (j < match->adjacency_total[edge.node]);
          )
          {
            iedge = match->clique_adjacency_list
              [match->clique.edge_total - 1][i];
            jedge = match->adjacency_list[edge.node][j];

            if (jedge.node < iedge.node)
              j++;

            else if (iedge.node < jedge.node)
              i++;

            else
            {
              match->clique_adjacency_list
                [match->clique.edge_total]
                [match->clique_adjacency_total[match->clique.edge_total]++] =
                (iedge.residual > jedge.residual) ? iedge : jedge;

              i++, j++;
            }
          }
        }

/*
*       If this node fails STAGE2 check, then for the sake of degeneracy
*       checking, record that this clique *could* have existed
*       11/96 te
*/
        else if (match->degeneracy_flag)
          match->degenerate[match->clique.edge_total] = TRUE;
      }
    }

/*
*   If we cannot expand the current clique, then process it.
*   10/96 te
*/
    else if (match->clique.edge_total > 0)
    {
/*
*     Compute residual of clique and check its validity
*     1/97 te
*/
      for (i = match->clique.residual = 0; i < match->clique.edge_total; i++)
        match->clique.residual =
          MAX (match->clique.residual, match->clique.edge[i].residual);

/*
*     This clique is sufficient for generating an orientation, if:
*     1. The clique is large enough, AND
*     2. The clique passes the degeneracy filter
*        a. The filter is not requested, OR
*        b. The filter indicates the cliqe is not degenerate
*     11/96 te
*/
      if
      (
        (match->clique.edge_total >= match->clique_size_min) &&
        (!match->degeneracy_flag ||
          !match->degenerate[match->clique.edge_total]) &&
        (match->clique.residual >
          ((match->cycle - 1) * match->distance_tolerance))
      )
        valid_clique = TRUE;

/*
*     Otherwise, erase the last node and decrement the clique size
*     11/96 te
*/
      else
      {
        match->clique.edge_total--;
        match->clique.edge[match->clique.edge_total].node = 0;
        match->clique.edge[match->clique.edge_total].residual = 0;
      }
    }

    else
      return EOF;
  }

/*
* Update match statistics
* 11/96 te
*/
  match->histogram[match->clique.edge_total - 1]++;

/*
* Reset ligand reflection state
* 1/97 te
*/
  if (match->clique.edge_total < 4)
    match->clique.reflect_flag = FALSE;

  return TRUE;
}


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

Routine to check if the ligand and receptor portions of a clique
have the same or opposite chirality.

Return values:
	TRUE:	The ligand must be reflected
	FALSE:	The ligand must not be reflected
7/95 te

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

int match_chirality (MATCH *match, int node)
{
  int i;
  static int atoms[4], sites[4];

  float calc_chirality (int *, XYZ **);
  
  for (i = 0; i < 3; i++)
  {
    atoms[i] = match->clique.edge[i].node % match->ligand_center.total.atoms;
    sites[i] = match->clique.edge[i].node / match->ligand_center.total.atoms;
  }

  atoms[3] = node % match->ligand_center.total.atoms;
  sites[3] = node / match->ligand_center.total.atoms;

  if (calc_chirality (atoms, match->ligand_vectors) *
    calc_chirality (sites, match->receptor_vectors) < 0.0)
    return TRUE;

  else
    return FALSE;
}

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

Routine to calculate the whether the fourth member is above or
below the plane defined by the first three members of a group.
7/95 te

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

float calc_chirality (int atoms[4], XYZ **vector)
{
  static XYZ normal;
  void vcrossv3 (XYZ, XYZ, XYZ);
  float vdotv3 (XYZ, XYZ);

  vcrossv3 (vector[atoms[0]][atoms[1]], vector[atoms[0]][atoms[2]], normal);
  return vdotv3 (vector[atoms[0]][atoms[3]], normal);
}



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

Routine to output match routine statistics
7/95 te

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

void output_match_info (MATCH *match)
{
  int i;

/*
* Output histogram of clique formation
* 6/95 te
*/
  if ((global.output_volume != 't') && (match->total > 0))
  {
    fprintf (global.outfile, "\n_______Matching_Histogram______\n");

    for (i = 0; i < match->clique_size_max; i++)
      if (match->histogram[i] > 0)
        fprintf (global.outfile, "%7d node%s: %9d\n",
          i + 1, (i == 0 ? " " : "s"),
          match->histogram[i]);

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

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

void extract_clique (MATCH *match, MOLECULE *molecule)
{
  int i, lig, rec;

  molecule->transform.refl_flag = match->clique.reflect_flag;

  for (i = 0; i < match->clique.edge_total; i++)
  {
    rec = match->clique.edge[i].node / match->ligand_center.total.atoms;

    copy_atom
      (&match->receptor_clique.atom[i],
      &match->receptor_site.atom[rec]);
    copy_coord
      (match->receptor_clique.coord[i],
      match->receptor_site.coord[rec]);
  }

  for (i = 0; i < match->clique.edge_total; i++)
  {
    lig = match->clique.edge[i].node % match->ligand_center.total.atoms;

    copy_atom
      (&match->ligand_clique.atom[i],
      &match->ligand_center.atom[lig]);
    copy_coord
      (match->ligand_clique.coord[i],
      match->ligand_center.coord[lig]);
  }

  match->ligand_clique.total.atoms =
    match->receptor_clique.total.atoms = match->clique.edge_total;
}


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

Routine to construct a random match
3/97 te

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

int get_random_match
(
  MATCH		*match,
  LABEL		*label,
  MOLECULE      *mol_conf,
  int           molecule_id,
  int           conformation_id,
  int           orientation_id
)
{
  int i, j;
  int unique;

  if (orientation_id == 0)
  {
    if ((molecule_id == 0) && (conformation_id == 0))
    {
      get_site (match, label);
      get_centers (match, label);

      match->clique.edge_max = match->clique_size_min;
      allocate_clique (&match->clique);
      match->clique.edge_total = match->clique_size_min;

      allocate_clique_atoms (match);
    }
    
    if (!match->centers_flag)
      get_ligand_centers (match, mol_conf);

    if (conformation_id == 0)
      match->node_total =
        match->ligand_center.total.atoms * match->receptor_site.total.atoms;
  }

/*
* Assign nodes to clique
* 3/97 te
*/
  for (i = 0; i < match->clique.edge_total;)
  {
    match->clique.edge[i].node = rand () % match->node_total;

    for (j = 0, unique = TRUE; j < i; j++)
      if ((match->clique.edge[i].node % match->ligand_center.total.atoms ==
        match->clique.edge[j].node % match->ligand_center.total.atoms) ||
        (match->clique.edge[i].node / match->ligand_center.total.atoms ==
        match->clique.edge[j].node / match->ligand_center.total.atoms))
        unique = FALSE;

    if (unique)
      i++;
  }

  return TRUE;
}


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

Routine to free random matching arrays
5/97 te

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

void free_random_matches (MATCH *match)
{
  free_clique (&match->clique);
  free_ligand_centers (match);
  free_molecule (&match->receptor_site);
  free_molecule (&match->ligand_clique);
  free_molecule (&match->receptor_clique);
}


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

Routine to manage chemical screening
3/96 te

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

int check_screen
(
  MATCH		*match,
  LABEL		*label,
  MOLECULE      *molecule,
  int           molecule_id
)
{
  static MOLECULE temporary = {0};

  match->chiral_flag = FALSE;

  if (label->chemical.screen.pharmaco_flag)
  {
    if (molecule_id == 0)
    {
      init_match_site (match, label);
      update_keys (&label->chemical, &match->receptor_site, 0);
    }

    if (label->chemical.screen.fold_flag)
      return
        check_pharmacophore
        (
          &label->chemical,
          match->distance_tolerance,
          &match->receptor_site,
          molecule
        );

    else
    {
      get_ligand_keys (match, molecule);
      init_match_ligand (match, label);
      initialize_adjacency (match, label);

      if (match_distance (match) != EOF)
        return TRUE;

      else
        return FALSE;
    }
  }

  else if (label->chemical.screen.similar_flag)
  {
    if (molecule_id == 0)
      init_similar_site (match, label);

    molecule->score.total =
      molecule->score.inter.total =
      check_dissimilarity (&label->chemical, &match->receptor_site, molecule);

    if (molecule->score.total <=
      label->chemical.screen.dissimilar_maximum)
      return TRUE;

    else
      return FALSE;
  }

  else if (label->chemical.screen.fold_flag)
  {
    if (!label->chemical.init_flag)
      get_chemical_labels (&label->chemical);

    fold_keys (&label->chemical, &temporary, molecule);
    copy_keys (molecule, &temporary);
    return TRUE;
  }

  else
    return FALSE;
}



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

Routine to free chemical screening arrays
5/97 te

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

void free_screen
(
  MATCH		*match,
  LABEL		*label
)
{
  if (label->chemical.screen.pharmaco_flag)
  {
    free_match_site (match, label);
/*
    if (label->chemical.screen.fold_flag)
      return
        check_pharmacophore
        (
          &label->chemical,
          match->distance_tolerance,
          &match->receptor_site,
          molecule
        );

    else
    {
      get_ligand_keys (match, molecule);
      init_match_ligand (match, label);
      initialize_adjacency (match, label);

      if (match_distance (match) != EOF)
        return TRUE;

      else
        return FALSE;
    }
*/
  }


  else if (label->chemical.screen.similar_flag)
  {
    free_molecule (&match->receptor_site);
    free_table (&label->chemical, &label->chemical.screen_table);
  }

/*
  else if (label->chemical.screen.fold_flag)
  {
    if (!label->chemical.init_flag)
      get_chemical_labels (&label->chemical);

    fold_keys (&label->chemical, &temporary, molecule);
    copy_keys (molecule, &temporary);
    return TRUE;
  }
*/
}


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

Routine to prepare for similarity search
3/96 te

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

void init_similar_site
(
  MATCH *match,
  LABEL *label
)
{
  FILE *file;

  get_chemical_labels (&label->chemical);

/*
* Read in reference molecule for similarity search
* 11/96 te
*/
  file = efopen (match->receptor_file_name, "r", global.outfile);

  if (read_molecule
    (&match->receptor_site, NULL, match->receptor_file_name, file, FALSE)
    != TRUE)
    exit (fprintf (global.outfile,
      "ERROR init_similar_site: Error reading in reference molecule.\n"));

  fclose (file);
}