score.c

最经典的分子对结软件

/*                                                                    */
/*                        Copyright UCSF, 1997                        */
/*                                                                    */

/*
Written by Todd Ewing
10/95
*/

#include "define.h"
#include "utility.h"
#include "mol.h"
#include "global.h"
#include "label.h"
#include "score.h"
#include "io_grid.h"
#include "io_receptor.h"
#include "vector.h"

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

Routine to read in receptor atoms and partition them into a
3-D grid for rapid evaluation of continuous scores.

2/97 te

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

void make_receptor_grid
(
  SCORE_GRID		*grid,
  SCORE_ENERGY		*energy,
  LABEL			*label
)
{
  FILE *receptor_file;
  XYZ span;
  int i, j, index, grid_coord[3];
  SLINT **grid_current = NULL;
  SLINT **grid_previous = NULL;

  grid->init_flag = TRUE;

/*
* Open and read receptor atom file
* 9/96 te
*/
  receptor_file = efopen (grid->receptor_file_name, "r", global.outfile);

  read_receptor
  (
    energy,
    label,
    &grid->receptor,
    grid->receptor_file_name,
    receptor_file,
    label->vdw.flag || label->chemical.flag,
    label->chemical.flag,
    label->vdw.flag
  );

/*
* Determine size of box enclosing receptor
* 9/96 te
*/
  for (j = 0; j < 3; j++)
  {
    grid->origin[j] = grid->receptor.coord[0][j];
    span[j] = grid->receptor.coord[0][j];
  }

  for (i = 1; i < grid->receptor.total.atoms; i++)
    for (j = 0; j < 3; j++)
    {
      if (grid->receptor.coord[i][j] < grid->origin[j])
        grid->origin[j] = grid->receptor.coord[i][j];

      if (grid->receptor.coord[i][j] > span[j])
        span[j] = grid->receptor.coord[i][j];
    }

  grid->size = 1;

  for (j = 0; j < 3; j++)
  {
    grid->span[j] = NINT ((span[j] - grid->origin[j]) / grid->spacing) + 1;
    grid->size *= grid->span[j];
  }

  ecalloc
  (
    (void **) &grid->atom,
    grid->size,
    sizeof (SLINT *),
    "receptor atom grid",
    global.outfile
  );

  ecalloc
  (
    (void **) &grid_current,
    grid->size,
    sizeof (SLINT *),
    "receptor atom ptr grid",
    global.outfile
  );

  ecalloc
  (
    (void **) &grid_previous,
    grid->size,
    sizeof (SLINT *),
    "receptor atom ptr grid",
    global.outfile
  );

  for (i = 0; i < grid->receptor.total.atoms; i++)
  {
    for (j = 0; j < 3; j++)
      grid_coord[j] = NINT ((grid->receptor.coord[i][j] - grid->origin[j])
        / grid->spacing);

    index =
      grid->span[0] * grid->span[1] * grid_coord[2] +
      grid->span[0] * grid_coord[1] +
      grid_coord[0];

    ecalloc
    (
      (void **) &grid_current[index],
      1,
      sizeof (SLINT),
      "next atom in receptor atom grid",
      global.outfile
    );

    grid_current[index]->value = i;

    if (grid->atom[index])
      grid_previous[index]->next = grid_current[index];

    else
      grid->atom[index] = grid_current[index];

    grid_previous[index] = grid_current[index];
    grid_current[index] = NULL;
  }

  efree ((void **) &grid_current);
  efree ((void **) &grid_previous);
}


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

Routine to free receptor grid

5/97 te

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

void free_receptor_grid (SCORE_GRID *grid)
{
  int i;
  SLINT *previous = NULL;
  SLINT *current = NULL;

  free_molecule (&grid->receptor);

  if (grid->atom)
  {
    for (i = 0; i < grid->size; i++)
    {
      for
      (
        previous = grid->atom[i];
        previous != NULL;
        previous = current
      )
      {
        current = previous->next;
        efree ((void **) &previous);
      }
    }

    efree ((void **) &grid->atom);
  }
}


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

Routine to compute the 1D array index given 3D coordinates.

Return values:
	POSITIVE INTEGER	grid index
	NEITHER (-1)		coordinates outside grid

2/97 te

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

int get_grid_index
(
  SCORE_GRID	*grid,
  XYZ		coord
)
{
  int grid_coord[3];

  if (get_grid_coordinate (grid, coord, grid_coord))
    return
      grid->span[0] * grid->span[1] * grid_coord[2] +
      grid->span[0] * grid_coord[1] +
      grid_coord[0];

  else
    return NEITHER;
}


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

Routine to compute the 3D integer grid coordinates given
3D real coordinates.

Return values:
	TRUE		coordinates inside grid
	FALSE		coordinates outside grid

2/97 te

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

int get_grid_coordinate
(
  SCORE_GRID	*grid,
  XYZ		coord,
  int		grid_coord[3]
)
{
  int i;
  int inside_flag = TRUE;

  for (i = 0; i < 3; i++)
  {
    grid_coord[i] = NINT ((coord[i] - grid->origin[i]) / grid->spacing);

    if
    (
      (grid_coord[i] < 0) ||
      (grid_coord[i] > grid->span[i] - 1)
    )
      inside_flag = FALSE;
  }

  return inside_flag;
}


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

Routine to identify all receptor atoms near the ligand atom
and compute a continuous intermolecular score.

2/97 te

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

void calc_inter_score_cont
(
  SCORE_GRID	*grid,
  void		*score,
  float		distance_cutoff,
  void		calc_inter_score
                  (SCORE_GRID *, void *, LABEL *, MOLECULE *,
                  int, int, SCORE_PART *),
  LABEL		*label,
  MOLECULE	*molecule,
  int		atom,
  SCORE_PART	*inter
)
{
  int i, j, k;
  int ilo, jlo, klo;
  int ihi, jhi, khi;
  int grid_cutoff;
  int grid_coord[3];
  int index;
  int rec_atom;
  SLINT *ptr;

  get_grid_coordinate (grid, molecule->coord[atom], grid_coord);
  grid_cutoff = (int) (distance_cutoff / grid->spacing) + 1,

  ilo = MAX (0, (grid_coord[0] - grid_cutoff));
  jlo = MAX (0, (grid_coord[1] - grid_cutoff));
  klo = MAX (0, (grid_coord[2] - grid_cutoff));
  ihi = MIN (grid->span[0], (grid_coord[0] + grid_cutoff + 1));
  jhi = MIN (grid->span[1], (grid_coord[1] + grid_cutoff + 1));
  khi = MIN (grid->span[2], (grid_coord[2] + grid_cutoff + 1));

  for (i = ilo; i < ihi; i++)
    for (j = jlo; j < jhi; j++)
      for (k = klo; k < khi; k++)
      {
        index =
          grid->span[0] * grid->span[1] * k +
          grid->span[0] * j + i;

        for (ptr = grid->atom[index]; ptr; ptr = ptr->next)
        {
          rec_atom = ptr->value;

          calc_inter_score
          (
            grid,
            score,
            label,
            molecule,
            atom,
            rec_atom,
            inter
          );
        }
      }
}


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

Routine to evaluate the number of intermolecular bumps between
a ligand and receptor

2/97 te

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

int check_bump
(
  SCORE_GRID *grid,
  SCORE_BUMP *bump,
  LABEL *label,
  MOLECULE *molecule
)
{
  int atom;
  int segment;
  int index;

  if (bump->grid == NULL)
    exit (fprintf (global.outfile, "ERROR check_bump: No bump grid loaded\n"));

  for (atom = molecule->score.bumpcount = 0;
    (molecule->score.bumpcount <= bump->maximum) && 
    (atom < molecule->total.atoms); atom++)
  {
    if
    (
      ((segment = molecule->atom[atom].segment_id) != NEITHER) &&
      (segment < molecule->total.segments) &&
      (molecule->segment[segment].active_flag != TRUE)
    )
      continue;

    if (molecule->atom[atom].heavy_flag == TRUE)
    {
      index = get_grid_index (grid, molecule->coord[atom]);

      if (index != NEITHER)
      {
        if (grid->version < 3.99)
        {
          if (bump->grid[index] == 'T')
            molecule->score.bumpcount++;
        }

        else
        {
          if (label->vdw.member[molecule->atom[atom].vdw_id].bump_id >=
            bump->grid[index])
            molecule->score.bumpcount++;
        }
      }
    }
  }

  return molecule->score.bumpcount;
}


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

Routine to compute the intermolecular contact score between
a ligand atom and the receptor.

2/97 te

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

void calc_inter_contact
(
  SCORE_GRID	*grid,
  SCORE_BUMP	*bump,
  SCORE_CONTACT	*contact,
  LABEL		*label,
  MOLECULE	*molecule,
  int		atom,
  SCORE_PART	*inter
)
{
  if (molecule->atom[atom].heavy_flag == TRUE)
  {
    if (grid->flag)
      calc_inter_contact_grid
        (grid, bump, contact, label, molecule, atom, inter);

    else
      calc_inter_score_cont
      (
        grid,
        (void *) contact,
        contact->distance,
        (void (*)()) calc_inter_contact_cont,
        label,
        molecule,
        atom,
        inter
      );
  }
}


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

Routine to compute the intermolecular contact score using a 
precomputed grid.

2/97 te

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

void calc_inter_contact_grid
(
  SCORE_GRID	*grid,
  SCORE_BUMP	*bump,
  SCORE_CONTACT	*contact,
  LABEL		*label,
  MOLECULE	*molecule,
  int		atom,
  SCORE_PART	*inter
)
{
  int index;

  index = get_grid_index (grid, molecule->coord[atom]);

  if (index != NEITHER)
  {
    if (label->vdw.member[molecule->atom[atom].vdw_id].bump_id >=
      bump->grid[index])
      inter->total += contact->clash_penalty;

    else
      inter->total += (float) contact->grid[index];
  }
}


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

Routine to compute the intermolecular contact score in a continuous
fashion given a ligand atom and a receptor atom.

2/97 te

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

void calc_inter_contact_cont
(
  SCORE_GRID	*grid,
  SCORE_CONTACT	*contact,
  LABEL		*label,
  MOLECULE	*molecule,
  int		atom,
  int		rec_atom,
  SCORE_PART	*inter
)
{
  float score;

  if (grid->receptor.atom[rec_atom].heavy_flag == TRUE)
  {
    calc_pairwise_contact