score_dock.c

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

void write_periph_structures
(
  DOCK		*dock,
  SCORE		*score,
  LIST		*list,
  MOLECULE	*mol_ref,
  MOLECULE	*mol_ori,
  int		type,
  int		anchor,
  int		layer,
  int		segment,
  int		shrink_flag
)
{
  FILE_NAME original_name;	/* File name of output file */
  FILE *original_file;		/* File pointer of output file */
  int number_written;		/* Number of structures written to file */

/*
* Store, then edit name of output file
* 2/98 te
*/
  strcpy (original_name, score->type[type].file_name);
  original_file = score->type[type].file;
  number_written = score->type[type].number_written;
  score->type[type].number_written = 0;

  sprintf (strrchr (score->type[type].file_name, '.') , "-%d-%02d-%d-%d%s",
    anchor + 1, layer + 1, segment + 1, shrink_flag + 1,
    strrchr (original_name, '.'));
  score->type[type].file =
    efopen (score->type[type].file_name, "w", global.outfile);

  write_topscorers
  (
    dock,
    score,
    list,
    mol_ref,
    mol_ori
  );

  score->type[type].number_written = number_written;
  efclose (&score->type[type].file);
  strcpy (score->type[type].file_name, original_name);
  score->type[type].file = original_file;
}



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

void minimize_ligand
(
  DOCK		*dock,
  LABEL		*label,
  SCORE		*score,
  MOLECULE	*mol_ref,
  MOLECULE	*mol_ori,
  MOLECULE	*mol_score,
  int		rigid_flag,
  int		layer_inner,
  int		layer_outer
)
{
  int i, j;			/* Iteration variable */
  int torsion;			/* Torsion id */
  int segment;			/* Segment id */
  int layer;			/* Layer id */
  int neighbor;			/* Segment neighbor id */
  int flex;			/* Flex label id */
  SIMPLEX simplex = {0};	/* Simplex data */
  int iteration_total;		/* Number of simplex iterations */
  int vertex_total = 0;		/* Total simplex vertices */
  float *vertex = NULL;		/* Simplex vertices */
  float delta_cycle;		/* Change in score during a cycle */
  float delta_total;		/* Total change in score */
  float distance;		/* Simplex vertex distance */

  if ((score->intra_flag) && (label->flex.minimize_flag == TRUE))
  {
    if
    (
      (layer_inner < 0) ||
      (layer_inner > layer_outer) ||
      (layer_outer >= mol_score->total.layers)
    )
      exit (fprintf (global.outfile,
        "ERROR minimize_ligand: molecule layer selection inappropriate\n"));

/*
*   Identify inner neighbors of minimizing segments
*   4/97 te
*/
    for (layer = layer_outer; layer >= layer_inner; layer--)
    {
      for (i = 0; i < mol_score->layer[layer].segment_total; i++)
      {
        segment = mol_score->layer[layer].segment[i];

        if (mol_score->segment[segment].min_flag)
        {
          if
          (
            ((torsion = mol_score->segment[segment].torsion_id) != NEITHER) &&
            ((flex = mol_score->torsion[torsion].flex_id) != NEITHER) &&
            (label->flex.member[flex].minimize_flag == TRUE)
          )
          {
            vertex_total++;
            mol_score->transform.tors_flag = TRUE;
          }

          if (layer > layer_inner)
          {
            for (j = 0; j < mol_score->segment[segment].neighbor_total; j++)
            {
              neighbor = mol_score->segment[segment].neighbor[j].id;

              if
              (
                (mol_score->segment[neighbor].layer_id < layer) &&
                (mol_score->segment[neighbor].layer_id >= layer_inner)
              )
                mol_score->segment[neighbor].min_flag = TRUE;
            }
          }
        }
      }
    }

/*
*   If no inner bonds can be minimized, then perform rigid-body minimization
*   4/97 te
*/
    if (vertex_total == 0)
      rigid_flag = TRUE;

/*
*   Identify outer neighbors of minimizing segments (not already identified)
*   4/97 te
*/
    for (layer = layer_inner + 1; layer <= layer_outer; layer++)
    {
      for (i = 0; i < mol_score->layer[layer].segment_total; i++)
      {
        segment = mol_score->layer[layer].segment[i];

        if (mol_score->segment[segment].min_flag != TRUE)
        {
          for (j = 0; j < mol_score->segment[segment].neighbor_total; j++)
          {
            neighbor = mol_score->segment[segment].neighbor[j].id;

            if
            (
              (mol_score->segment[neighbor].layer_id < layer) &&
              (mol_score->segment[neighbor].layer_id >= layer_inner) &&
              (mol_score->segment[neighbor].min_flag == TRUE)
            )
            {
              if
              (
                ((torsion = mol_score->segment[segment].torsion_id)
                  != NEITHER) &&
                ((flex = mol_score->torsion[torsion].flex_id) != NEITHER) &&
                (label->flex.member[flex].minimize_flag == TRUE)
              )
              {
                vertex_total++;
                mol_score->transform.tors_flag = TRUE;
              }

              mol_score->segment[segment].min_flag = TRUE;
              break;
            }
          }
        }
      }
    }
  }

  if ((score->inter_flag) && (rigid_flag == TRUE))
  {
    simplex.rigid_flag = TRUE;
    vertex_total += 6;
    mol_score->transform.trans_flag =
      mol_score->transform.rot_flag = TRUE;
  }

  if (vertex_total > 0)
  {
    ecalloc
    (
      (void **) &vertex,
      vertex_total,
      sizeof (float),
      "simplex vertex array",
      global.outfile
    );

    simplex.dock = dock;
    simplex.score = score;
    simplex.label = label;
    simplex.mol_ref = mol_ref;
    simplex.mol_ori = mol_ori;
    simplex.mol_score = mol_score;
    simplex.layer_inner = layer_inner;
    simplex.layer_outer = layer_outer;

    copy_molecule (&simplex.mol_min, mol_score);
    copy_molecule (&simplex.mol_best, mol_score);
    simplex.mol_best.score.total = INITIAL_SCORE;
    delta_total = 0;

    for
    (
      simplex.cycle = 0, distance = FLT_MAX;
      (simplex.cycle < score->minimize.cycle) &&
        (distance > score->minimize.cycle_converge) &&
        ((simplex.cycle == 0) ||
          (simplex.mol_best.score.total <
          score->type[mol_score->score.type].termination));
      simplex.cycle++
    )
    {
      simplex_optimize
      (
        (void *) &simplex,
        vertex,
        vertex_total,
        score->type[mol_score->score.type].convergence,
        &iteration_total,
        score->minimize.iteration,
        simplex_score,
        &delta_cycle
      );

      copy_molecule (mol_score, &simplex.mol_best);

/*
*     Compute distance that vertex has moved (and update initial vertex)
*     1/97 te
*/
      for (i = distance = 0; i < vertex_total; i++)
        distance += SQR (vertex[i]);

      distance = sqrt (distance) / (float) (simplex.cycle + 1);
      memset (vertex, 0, vertex_total * sizeof (float));

/*
*     Update minimizer statistics
*     1/97 te
*/
      delta_total += delta_cycle;

      score->minimize.iteration_total += iteration_total;
      score->minimize.iteration_max =
        MAX (score->minimize.iteration_max, iteration_total);
      score->minimize.iteration_min =
        MIN (score->minimize.iteration_min, iteration_total);
    }

    score->minimize.call_sub_total++;

    score->minimize.vertex_total += vertex_total;
    score->minimize.vertex_max =
      MAX (score->minimize.vertex_max, vertex_total);
    score->minimize.vertex_min =
      MIN (score->minimize.vertex_min, vertex_total);

    score->minimize.cycle_total += simplex.cycle;
    score->minimize.cycle_max =
      MAX (score->minimize.cycle_max, simplex.cycle);
    score->minimize.cycle_min =
      MIN (score->minimize.cycle_min, simplex.cycle);

    score->minimize.delta_total += delta_total;
    score->minimize.delta_max =
      MAX (score->minimize.delta_max, delta_total);
    score->minimize.delta_min =
      MIN (score->minimize.delta_min, delta_total);

    free_molecule (&simplex.mol_min);
    free_molecule (&simplex.mol_best);
    efree ((void **) &vertex);
  }

  else
    calc_score
    (
      label,
      score,
      mol_score,
      layer_outer
    );

/*
* Turn off all minimize flags
* 3/97 te
*/
  if (label->flex.minimize_flag == TRUE)
    for (layer = layer_inner; layer <= layer_outer; layer++)
      for (i = 0; i < mol_score->layer[layer].segment_total; i++)
      {
        segment = mol_score->layer[layer].segment[i];
        mol_score->segment[segment].min_flag = FALSE;
      }
}

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

float simplex_score (void *simplex_input, float *vertex)
{
  int i;
  int flex;			/* Flex label id */
  int torsion;			/* Torsion id */
  int segment;			/* Segment id */
  int layer;			/* Layer id */
  int intra_current_flag;	/* Flag for whether score up-to-date */
  int vertex_count = 0;
  SIMPLEX *simplex;
  float score;

  simplex = (SIMPLEX *) simplex_input;

/*
* Extract transformation variables from simplex array if zeroth layer
* 12/96 te
*/
  if ((simplex->score->inter_flag) && (simplex->rigid_flag == TRUE))
  {
    for (i = 0; i < 3; i++)
    {
      simplex->mol_min.transform.translate[i] =
        simplex->mol_score->transform.translate[i] +
        vertex[i] * simplex->score->minimize.translation /
        (float) (simplex->cycle + 1);

      simplex->mol_min.transform.rotate[i] =
        simplex->mol_score->transform.rotate[i] +
        vertex[i + 3] * simplex->score->minimize.rotation /
        (float) (simplex->cycle + 1);
    }

    transform_molecule (&simplex->mol_min, simplex->mol_ref);

    vertex_count = 6;

/*
*   Flag all segment intermolecular scores as out-of-date
*   3/97 te
*/
    for (segment = 0; segment < simplex->mol_min.total.segments; segment++)
      simplex->mol_min.segment[segment].score.inter.current_flag = FALSE;
  }

  else
    vertex_count = 0;

/*
* Extract torsion variables from simplex array
* 3/97 te
*/
  if ((simplex->score->intra_flag) &&
    (simplex->label->flex.minimize_flag == TRUE))
  {
    intra_current_flag = TRUE;

    for (layer = simplex->layer_inner; layer <= simplex->layer_outer; layer++)
    {
      for (i = 0; i < simplex->mol_min.layer[layer].segment_total; i++)
      {
        segment = simplex->mol_min.layer[layer].segment[i];

        if ((simplex->mol_min.segment[segment].min_flag == TRUE) &&
          ((torsion = simplex->mol_min.segment[segment].torsion_id) != NEITHER))
        {
          if
          (
            ((flex = simplex->mol_min.torsion[torsion].flex_id) != NEITHER) &&
            (simplex->label->flex.member[flex].minimize_flag == TRUE)
          )
          {
            simplex->mol_min.torsion[torsion].target_angle =
              simplex->mol_score->torsion[torsion].target_angle +
              vertex[vertex_count++] *
                simplex->score->minimize.torsion * PI / 180.0 /
                (float) (simplex->cycle + 1) /
                (float) (simplex->layer_outer - layer + 1);

            torsion_transform (&simplex->mol_min, torsion);
          }

          intra_current_flag = FALSE;
          simplex->mol_min.segment[segment].score.inter.current_flag = FALSE;
        }
      }

      if (intra_current_flag == FALSE)
        for (i = 0; i < simplex->mol_min.layer[layer].segment_total; i++)
        {
          segment = simplex->mol_min.layer[layer].segment[i];
          simplex->mol_min.segment[segment].score.intra.current_flag = FALSE;
        }
    }
  }

/*
* Evaluate the score of this ligand position
* 12/96 te
*/
  score =
    calc_score
    (
      simplex->label,
      simplex->score,
      &simplex->mol_min,
      simplex->layer_outer
    );

  if (score < simplex->mol_best.score.total)
    copy_molecule (&simplex->mol_best, &simplex->mol_min);

  return score;
}


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

float calc_score
(
  LABEL		*label,
  SCORE		*score,
  MOLECULE	*molecule,
  int		layer_outer
) 
{
  int li, lj;			/* Layer id */
  int sli, slj;			/* Segment id in layer list */
  int si, sj;			/* Segment id */
  int asi, asj;			/* Atom id in segment list */
  int ai, aj;			/* Atom id */

/*
* Update near flag array which speeds up intramolecular score calculation
* 3/97 te
*/
  if (score->intra_flag)
    initialize_near (&score->near, molecule);

/*
* Initialize score data
* 1/97 te
*/
  reset_score_parts (&molecule->score.intra);
  reset_score_parts (&molecule->score.inter);
  molecule->score.total = 0;

/*
* Loop through all layers of interest
* 1/97 te
*/
  for (li = layer_outer; li >= 0; li--)
  {
/*
*   Loop through layer segments
*   3/97 te
*/
    for (sli = molecule->layer[li].segment_total - 1; sli >= 0; sli--)
    {
      si = molecule->layer[li].segment[sli];
      if (!molecule->segment[si].active_flag) continue;

/*
*     Calculate intermolecular score
*     3/97 te
*/
      if (score->inter_flag)
      {
        if (!molecule->segment[si].score.inter.current_flag ||
          (molecule->segment[si].score.type != molecule->score.type))
        {
          reset_score_parts (&molecule->segment[si].score.inter);
          molecule->segment[si].score.inter.current_flag = TRUE;

/*
*         Loop through segment atoms
*         3/97 te
*/
          for (asi = molecule->segment[si].atom_total - 1; asi >= 0; asi--)
          {
            ai = molecule->segment[si].atom[asi];

            calc_inter_score
            (
              label,
              score,
              molecule,
              ai,
              &molecule->segment[si].score.inter
            );
          }
        }

        sum_score
        (
          molecule->score.type,
          &molecule->score.inter,
          &molecule->segment[si].score.inter
        );
      }

/*