score.c

最经典的分子对结软件

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

void calc_inter_chemical
(
  SCORE_GRID	*grid,
  SCORE_ENERGY	*energy,
  SCORE_CHEMICAL *chemical,
  LABEL		*label,
  MOLECULE	*molecule,
  int		atom,
  SCORE_PART	*inter
)
{
  if (!energy->vdw_init_flag)
    initialize_vdw_energy (energy, &label->vdw);

  if (label->chemical.score_table == NULL)
    get_table
    (
      &label->chemical,
      label->chemical.score_file_name,
      &label->chemical.score_table
    );

  if (label->vdw.member[molecule->atom[atom].vdw_id].well_depth != 0.0)
  {
    if (grid->flag)
      calc_inter_chemical_grid
        (grid, energy, chemical, label, molecule, atom, inter);

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


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

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

2/97 te

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

void calc_inter_chemical_grid
(
  SCORE_GRID	*grid,
  SCORE_ENERGY	*energy,
  SCORE_CHEMICAL *chemical,
  LABEL		*label,
  MOLECULE	*molecule,
  int		atom,
  SCORE_PART	*inter
)
{
  int i;
  int index[8];
  int out_flag;
  XYZ cube_coord;
  float vdw, electro;

  void get_index_ (XYZ, XYZ, float *, int *, int *, int *, float *);
  float get_value_ (float *, XYZ, int *);

  if (label->vdw.member[molecule->atom[atom].vdw_id].well_depth != 0.0)
  {
    get_index_
    (
      molecule->coord[atom],
      grid->origin,
      &grid->spacing,
      grid->span,
      &out_flag,
      index,
      cube_coord
    );

    if (out_flag == 0)
    {
      vdw =
        energy->vdwA[molecule->atom[atom].vdw_id] *
        get_value_(energy->avdw, cube_coord, index);

      for (i = 0; i < label->chemical.total; i++)
      {
        if (chemical->grid[i] != NULL)
          vdw -=
            label->chemical.score_table[molecule->atom[atom].chem_id][i] *
            energy->vdwB[molecule->atom[atom].vdw_id] *
            get_value_(chemical->grid[i], cube_coord, index);
      }

      electro =
        molecule->atom[atom].charge *
        get_value_(energy->es, cube_coord, index);
    }

    else
    {
      vdw = 10.0;
      electro = 10.0;
    }

    inter->vdw += vdw * energy->scale_vdw;
    inter->electro += electro * energy->scale_electro;
    inter->total = inter->vdw + inter->electro;
  }
}


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

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

2/97 te

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

void calc_inter_chemical_cont
(
  SCORE_GRID	*grid,
  SCORE_ENERGY	*energy,
  LABEL		*label,
  MOLECULE	*molecule,
  int		atom,
  int		rec_atom,
  SCORE_PART	*inter
)
{
  float vdwA, vdwB, electro;

  if (label->vdw.member[grid->receptor.atom[rec_atom].vdw_id].well_depth != 0.0)
  {
    calc_pairwise_energy
    (
      energy,
      label,
      molecule,
      &grid->receptor,
      atom,
      rec_atom,
      &vdwA,
      &vdwB,
      &electro
    );

    inter->vdw += (vdwA - vdwB *
      label->chemical.score_table
        [molecule->atom[atom].chem_id]
        [grid->receptor.atom[rec_atom].chem_id]) * energy->scale_vdw;

    inter->electro += electro * energy->scale_electro;
    inter->total = inter->vdw + inter->electro;
  }
}


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

Routine to compute the intramolecular chemical score
between two ligand atoms.

2/97 te

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

void calc_intra_chemical
(
  SCORE_ENERGY	*energy,
  LABEL		*label,
  MOLECULE	*molecule,
  int		atomi,
  int		atomj,
  SCORE_PART	*intra
)
{
  float vdwA, vdwB, electro;

  calc_pairwise_energy
  (
    energy,
    label,
    molecule,
    molecule,
    atomi, atomj,
    &vdwA,
    &vdwB,
    &electro
  );

  intra->vdw += (vdwA - vdwB *
    label->chemical.score_table
      [molecule->atom[atomi].chem_id]
      [molecule->atom[atomj].chem_id]) * energy->scale_vdw;

  intra->electro += electro * energy->scale_electro;
  intra->total = intra->vdw + intra->electro;
}


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

Routine to add chemical score components.

2/97 te

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

void sum_chemical
(
  SCORE_PART	*sum,
  SCORE_PART	*increment
)
{
  sum->vdw += increment->vdw;
  sum->electro += increment->electro;
  sum->total = sum->vdw + sum->electro;
}


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

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

2/97 te

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

void calc_inter_rmsd
(
  SCORE_GRID	*grid,
  MOLECULE	*molecule,
  int		atom,
  SCORE_PART	*inter
)
{
/*
* Check if input is alright
* 1/97 te
*/
  if (molecule->total.atoms != grid->receptor.total.atoms)
    exit (fprintf (global.outfile,
      "ERROR calc_rmsd_score: molecules have different number of atoms\n"));

  if (molecule->atom[atom].heavy_flag == TRUE)
  {
    inter->vdw += 1;
    inter->electro +=
      square_distance (molecule->coord[atom], grid->receptor.coord[atom]);
  }
}


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

Routine to add rmsd score components.

2/97 te

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

void sum_rmsd
(
  SCORE_PART	*sum,
  SCORE_PART	*increment
)
{
  sum->vdw += increment->vdw;
  sum->electro += increment->electro;

  if (sum->vdw > 0)
    sum->total = sqrt (sum->electro / sum->vdw);

  else
    sum->total = 0;
}


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

Routine to compute the rms deviation between two configurations 
of the same molecule.

2/97 te

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

float calc_rmsd
(
  MOLECULE	*mol_ori,
  MOLECULE	*mol_ref
)
{
  int atom;
  int layer_ref;
  int layer_ori;
  int heavy_total = 0;
  float rmsd = 0;

  if (mol_ori->total.atoms != mol_ref->total.atoms)
    exit (fprintf (global.outfile,
      "ERROR calc_rmsd: molecules have different number of atoms\n"));

/*
* If both molecules have an anchor layer, but the anchors are
* different, then set the rmsd arbitrarily high
* 2/97 te
  if ((mol_ori->total.layers > 1) &&
    (mol_ref->total.layers > 1))
  {
    if (mol_ori->layer[0].segment_total != mol_ref->layer[0].segment_total)
      return FLT_MAX;

    for (segment = 0; segment < mol_ref->layer[0].segment_total; segment++)
      if (mol_ori->layer[0].segment[segment] !=
        mol_ref->layer[0].segment[segment])
        return FLT_MAX;
  }
*/

/*
* Loop through atoms
* 2/97 te
*/
  for (atom = 0; atom < mol_ori->total.atoms; atom++)
  {
    if (mol_ori->atom[atom].heavy_flag != TRUE)
      continue;

    if ((mol_ori->total.segments > 0) &&
      (mol_ori->segment[mol_ori->atom[atom].segment_id].active_flag == FALSE))
      continue;
    
    if ((mol_ref->total.layers > 1) && (mol_ori->total.layers > 1))
    {
      layer_ref = mol_ref->segment[mol_ref->atom[atom].segment_id].layer_id;
      layer_ori = mol_ori->segment[mol_ori->atom[atom].segment_id].layer_id;

      if
      (
        ((layer_ref == NEITHER) || (layer_ori == NEITHER)) &&
        (layer_ref != layer_ori)
      )
        return FLT_MAX;
    }

    rmsd += square_distance (mol_ori->coord[atom], mol_ref->coord[atom]);
    heavy_total++;
  }

  if (heavy_total > 0)
    rmsd = sqrt (rmsd / (float) heavy_total);

  else
    rmsd = 0;

  return rmsd;
}


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

Routine to compute the layer-weighted rms deviation between two configurations 
of the same molecule.

3/97 te

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

float calc_layer_rmsd
(
  MOLECULE	*mol_ori,
  MOLECULE	*mol_ref
)
{
  int atom;
  int layer_ref;
  int layer_ori;
  int heavy_total = 0;
  float rmsd = 0;

  if (mol_ori->total.atoms != mol_ref->total.atoms)
    exit (fprintf (global.outfile,
      "ERROR calc_rmsd: molecules have different number of atoms\n"));

/*
* Loop through atoms
* 3/97 te
*/
  for (atom = 0; atom < mol_ori->total.atoms; atom++)
  {
/*
    if (mol_ori->atom[atom].heavy_flag != TRUE)
      continue;
*/

    if ((mol_ori->total.segments > 0) &&
      (mol_ori->segment[mol_ori->atom[atom].segment_id].active_flag == FALSE))
      continue;
    
    layer_ref = mol_ref->segment[mol_ref->atom[atom].segment_id].layer_id;
    layer_ori = mol_ori->segment[mol_ori->atom[atom].segment_id].layer_id;

    if
    (
      ((layer_ref == NEITHER) || (layer_ori == NEITHER)) &&
      (layer_ref != layer_ori)
    )
      return FLT_MAX;

    rmsd +=
      square_distance (mol_ori->coord[atom], mol_ref->coord[atom]) *
      (layer_ori + 1);
    heavy_total += layer_ori + 1;
  }

  if (heavy_total > 0)
    rmsd = sqrt (rmsd / (float) heavy_total);

  else
    rmsd = 0;

  return rmsd;
}


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

Routine to compute the rms deviation of a segment from two configurations 
of the same molecule.

3/97 te

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

float calc_segment_rmsd
(
  MOLECULE	*mol_ori,
  MOLECULE	*mol_ref,
  int		segment
)
{
  int atom;
  int atom_id;
  int atom_total = 0;
  float rmsd = 0;

  if (mol_ori->total.atoms != mol_ref->total.atoms)
    exit (fprintf (global.outfile,
      "ERROR calc_segment_rmsd: molecules have different number of atoms\n"));

/*
* Loop through atoms
* 2/97 te
*/
  for (atom_id = 0; atom_id < mol_ori->segment[segment].atom_total; atom_id++)
  {
    atom = mol_ori->segment[segment].atom[atom_id];

    if ((atom < 0) || (atom >= mol_ori->total.atoms))
      exit (fprintf (global.outfile,
        "ERROR calc_segment_rmsd: segment contains bad atom\n"));

    rmsd += square_distance (mol_ori->coord[atom], mol_ref->coord[atom]);
    atom_total++;                                                       
  }                                                                      

  if (atom_total > 0)
    rmsd = sqrt (rmsd / (float) atom_total);

  else                                       
    rmsd = 0;

  return rmsd;
}