match.c

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

Routine to set up arrays used during matching.
7/95 te

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

int initialize_adjacency
(
  MATCH		*match,
  LABEL		*label
)
{
  reset_match (match);

  if (label->chemical.screen.process_flag)
    compute_chemical_adjacency (match, label);

  else
    compute_adjacency (match);

  match->clique_adjacency_total[0] = match->node_total;

  return TRUE;
}


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

Routine to get site points
11/96 te

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

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

/*
* Read in receptor site points
* 3/96 te
*/
  file = efopen (match->receptor_file_name, "r", global.outfile);

  if (read_molecule
    (&match->receptor_site, NULL, match->receptor_file_name, file, TRUE)
    != TRUE)
    exit (fprintf (global.outfile,
      "ERROR get_site: Unable to read receptor site points from %s\n",
      match->receptor_file_name));

  if (match->clique_size_min > match->receptor_site.total.atoms)
    exit (fprintf (global.outfile,
      "ERROR get_site: clique_size_min > number of site points.\n"));

  if (prepare_molecule
    (
      &match->receptor_site,
      match->receptor_file_name,
      file,
      label,
      'a',
      FALSE,
      match->chemical_flag,
      label->chemical.screen.flag
    ) != TRUE)
    exit (fprintf (global.outfile,
      "ERROR get_site: Unable to prepare receptor site points.\n"));

  fclose (file);
}


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

Routine to get (or just allocate) ligand centers
3/96 te

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

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

/*
* If requested, read in ligand centers
* 3/96 te
*/
  if (match->centers_flag)
  {
    file = efopen (match->ligand_file_name, "r", global.outfile);

    if (read_molecule
      (
        &match->ligand_center,
        NULL,
        match->ligand_file_name,
        file,
        TRUE
      ) != TRUE)
      exit (fprintf (global.outfile,
        "ERROR get_centers: Unable to read ligand site points from %s\n",
        match->ligand_file_name));

    if (prepare_molecule
      (
        &match->ligand_center,
        match->ligand_file_name,
        file,
        label,
        'a',
        FALSE,
        match->chemical_flag,
        FALSE
      ) != TRUE)
      exit (fprintf (global.outfile,
        "ERROR get_centers: Unable to prepare ligand site points.\n"));

    center_of_mass
    (
      match->ligand_center.coord,
      match->ligand_center.total.atoms,
      match->ligand_center.transform.com
    );

    fclose (file);

    calculate_distances
    (
      &match->ligand_center,
      &match->ligand_distances,
      &match->ligand_distance_size
    );

    calculate_vectors
    (
      &match->ligand_center,
      &match->ligand_vectors,
      &match->ligand_vector_size
    );
  }
}


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

Routine to fix the critical cluster identifiers of the site points
so that the zeroth cluster (unclustered site points) becomes the
last cluster.
7/95 te

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

void id_site_points (MATCH *match)
{
  int i;

  if (match->critical_flag)
  {
    match->cluster_total = match->receptor_site.total.substs;

/*
*   Make sure the largest cluster is not greater than the largest
*   allowed clique size
*   7/95 te
*/
    if (match->cluster_total > match->clique_size_max)
      exit (fprintf (global.outfile,
        "ERROR id_site_points: "
        "Critical cluster identifiers cannot exceed nodes_maximum,\n"
        "      or if a zeroth cluster is present, nodes_maximum - 1.\n"));
  }

/*
* If critical cluster filtering is not to be performed, then
* set all cluster identifiers to zero.
* 7/95 te
*/
  else
  {
    for (i = 0; i < match->receptor_site.total.atoms; i++)
      match->receptor_site.atom[i].subst_id = 0;

    match->cluster_total = 0;
  }
}


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

Routine to reshuffle order of the site points to get the members of
critical clusters first
7/95 te

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

void order_site_points (MATCH *match)
{
  int i, j;
  MOLECULE temporary = {0};
  int compare_site_points ();

/*
* Allocate space for the cluster lists
* 11/96 te
*/
  ecalloc
  (
    (void **) &match->cluster,
    match->cluster_total,
    sizeof (int *),
    "cluster list",
    global.outfile
  );

  ecalloc
  (
    (void **) &match->cluster_size,
    match->cluster_total,
    sizeof (int),
    "cluster list",
    global.outfile
  );

/*
* Determine the size of each cluster
* 11/96 te
*/
  for (i = 0; i < match->receptor_site.total.atoms; i++)
    match->cluster_size[match->receptor_site.atom[i].subst_id]++;
  
  for (i = 0; i < match->cluster_total; i++)
    ecalloc
    (
      (void **) &match->cluster[i],
      match->cluster_size[i],
      sizeof (int),
      "cluster list",
      global.outfile
    );

/*
* Allocate memory for space to reorder site points
* 7/95 te
*/
  temporary.max.atoms = match->receptor_site.total.atoms;
  vstrcpy (&temporary.info.name, "temporary");
  allocate_molecule (&temporary);

/*
* Copy over site points in critical clusters
* 7/95 te
*/
  for (i = temporary.total.atoms = 0; i < match->cluster_total; i++)
    for (j = match->cluster_size[i] = 0;
      j < match->receptor_site.total.atoms; j++)
      if (match->receptor_site.atom[j].subst_id == i)
      {
        copy_atom
          (&temporary.atom[temporary.total.atoms],
          &match->receptor_site.atom[j]);
        copy_coord
          (temporary.coord[temporary.total.atoms],
          match->receptor_site.coord[j]);

        match->cluster[i][match->cluster_size[i]++] = temporary.total.atoms++;
      }

/*
* Check to make the correct number of atoms were copied
* 7/95 te
*/
  if (match->receptor_site.total.atoms != temporary.total.atoms)
    exit (fprintf (global.outfile,
      "ERROR order_site_points: "
      "Critical cluster identifiers must be in the range "
      "[0, nodes_maximum].\n"));

/*
* Copy the reordered site points back into molecular structure
* 7/95 te
*/
  copy_atoms (&match->receptor_site, &temporary);
  free_molecule (&temporary);

/*
* Check to make sure the first cluster is identified as number 1
* 7/95 te
*/
  if (match->receptor_site.atom[0].subst_id != 0)
    exit (fprintf (global.outfile,
      "ERROR order_site_points: "
      "Critical cluster identifiers must begin at 1.\n"));

/*
* Check to make sure the clusters are numbered sequentially with no gaps
* 7/95 te
*/
  for (i = 1; i < match->receptor_site.total.atoms; i++)
    if ((match->receptor_site.atom[i].subst_id !=
      match->receptor_site.atom[i - 1].subst_id) &&
      (match->receptor_site.atom[i].subst_id !=
      match->receptor_site.atom[i - 1].subst_id + 1))
      exit (fprintf (global.outfile,
        "ERROR order_site_points: "
        "Critical cluster identifiers cannot have gaps.\n"));

/*
* Print out critical clusters
* 7/95 te
*/
  fprintf (global.outfile, "____Critical_Clusters____\n");

  fprintf (global.outfile, "%5s %5s %5s\n", "clust", "name", "type");
  for (i = 0; (i < match->receptor_site.total.atoms) &&
    (match->receptor_site.atom[i].subst_id <= match->cluster_total); i++)
    fprintf (global.outfile, "%5d %5s %5s\n",
      match->receptor_site.atom[i].subst_id + 1,
      match->receptor_site.atom[i].name, match->receptor_site.atom[i].type);

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


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

Routine to free critical clusters
5/97 te

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

void free_site_points (MATCH *match)
{
  int i;

/*
* Free space for the cluster lists
* 11/96 te
*/
  for (i = 0; i < match->cluster_total; i++)
    efree ((void **) &match->cluster[i]);

  efree ((void **) &match->cluster);
  efree ((void **) &match->cluster_size);
}



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

Routine to construct a vector matrix: a N by N matrix in which each
element is a difference vector between the positions of each of N atoms
or site points.  This matrix is constructed to speed up the clique
chirality checking algorithms.
7/95 te

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

void calculate_vectors (MOLECULE *molecule, XYZ ***matrix, int *size)
{
  int i, j, k;

/*
* Either allocate space for difference vectors, or reset space
* 3/97 te
*/
  if (*size < molecule->total.atoms)
  {
    free_vectors (matrix, size);

    *size = molecule->total.atoms;

    ecalloc
    (
      (void **) matrix,
      *size,
      sizeof (XYZ *),
      "vectors matrix",
      global.outfile
    );

    for (i = 0; i < *size; i++)
      ecalloc
      (
        (void **) &(*matrix)[i],
        *size,
        sizeof (XYZ),
        "vectors matrix",
        global.outfile
      );
  }

  else
    for (i = 0; i < *size; i++)
      memset ((*matrix)[i], 0, *size * sizeof (XYZ));

/*
* Loop through all pairs of atoms/site points
* 3/97 te
*/
  for (i = 0; i < molecule->total.atoms; i++)
    for (j = 0; j < molecule->total.atoms; j++)
      for (k = 0; k < 3; k++)
        (*matrix)[i][j][k] =
          molecule->coord[j][k] -
          molecule->coord[i][k];
}


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

void free_vectors (XYZ ***matrix, int *size)
{
  int i;

  for (i = 0; i < *size; i++)
    efree ((void **) &(*matrix)[i]);

  efree ((void **) matrix);
  *size = 0;
}


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

Routine to construct the receptor distance bins used in match_driver.
1/97 te

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

void make_distance_bins (MATCH *match)
{
  int i, j, k;
  SLINT2 **current = NULL;
  SLINT2 **previous = NULL;

/*
* Allocate space for receptor distance bins
* 1/97 te
*/
  match->bin_width = 0.5;
  match->bin_total = NINT (match->bin_length / match->bin_width) + 1;

  ecalloc
  (
    (void **) &match->bin,
    match->bin_total,
    sizeof (SLINT2 *),
    "receptor distance bins",
    global.outfile
  );

  ecalloc
  (
    (void **) &previous,
    match->bin_total,
    sizeof (SLINT2 *),
    "receptor distance bins",
    global.outfile
  );

  ecalloc
  (
    (void **) &current,
    match->bin_total,
    sizeof (SLINT2 *),
    "receptor distance bins",
    global.outfile
  );

/*
* Divide receptor distances into bins
* 1/97 te
*/
  for (i = 0; i < match->receptor_site.total.atoms; i++)
    for (j = i + 1; j < match->receptor_site.total.atoms; j++)
    {
      k = NINT (match->receptor_distances[i][j] / match->bin_width);

      ecalloc
      (
        (void **) &current[k],
        1,
        sizeof (SLINT2),
        "next bin occupant",
        global.outfile
      );

      current[k]->i = i;
      current[k]->j = j;

      if (match->bin[k] != NULL)
        previous[k]->next = current[k];

      else
        match->bin[k] = current[k];

      previous[k] = current[k];
      current[k] = NULL;
    }

  efree ((void **) &current);
  efree ((void **) &previous);
}


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

Routine to free the receptor distance bins used in match_driver.
5/97 te

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

void free_distance_bins (MATCH *match)
{
  int i;
  SLINT2 *prev = NULL;
  SLINT2 *next = NULL;

  for (i = 0; i < match->bin_total; i++)
    for (prev = match->bin[i]; prev != NULL; prev = next)
    {
      next = prev->next;
      efree ((void **) &prev);
    }

  match->bin_total = 0;
  efree ((void **) &match->bin);
}


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

Routine to extract ligand centers from a ligand
3/96 te

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

int get_ligand_centers
(
  MATCH		*match,
  MOLECULE	*molecule
)
{
  int atom;
  int segment;
  int layer;

  match->ligand_center.total.atoms = molecule->total.atoms;
  reallocate_atoms (&match->ligand_center);

/*
* Allocate space for the ligand atom key
* 3/97 te
*/
  if (match->ligand_key_total < match->ligand_center.total.atoms)
  {
    efree ((void **) &match->ligand_key);
    match->ligand_key_total = match->ligand_center.total.atoms;

    ecalloc
    (
      (void **) &match->ligand_key,
      match->ligand_key_total,
      sizeof (int),
      "match ligand key",
      global.outfile
    );
  }

/*
* Copy the anchor heavy atoms over
* 3/96 te
*/
  for (atom = match->ligand_center.total.atoms = 0;
    atom < molecule->total.atoms; atom++)
  {
    if
    (
      ((segment = molecule->atom[atom].segment_id) != NEITHER) &&
      (segment < molecule->total.segments) &&
      ((layer = molecule->segment[segment].layer_id) != NEITHER) &&
      (layer < molecule->total.layers) &&
      (layer != 0)
    )
      continue;

    if (molecule->atom[atom].heavy_flag != TRUE)
      continue;

/*
    fprintf (global.outfile, "Lig cent %s, seg %d, layer %d\n",
      molecule->atom[atom].name, segment, layer);

    if (match->ligand_center.total.atoms >= match->ligand_center.max.atoms)
      exit (fprintf (global.outfile,
        "ERROR get_ligand_centers: Too many centers in molecule\n"));
*/

    copy_atom
      (&match->ligand_center.atom[match->ligand_center.total.atoms],
      &molecule->atom[atom]);
    copy_coord
      (match->ligand_center.coord[match->ligand_center.total.atoms],
      molecule->coord[atom]);

    match->ligand_key[match->ligand_center.total.atoms] = atom;
    match->ligand_center.total.atoms++;
  }