label_flex.c

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

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

/*
Written by Todd Ewing
12/96
*/
#include "define.h"
#include "utility.h"
#include "mol.h"
#include "global.h"
#include "search.h"
#include "label_node.h"
#include "label_vdw.h"
#include "label_flex.h"
#include "transform.h"


/* ==================================================================== */
int get_flex_labels (LABEL_FLEX *label_flex)
{
  int i, definition_total, definition_count;
  STRING100 line;
  FILE *flex_file;

  label_flex->init_flag = TRUE;

  flex_file = efopen (label_flex->file_name, "r", global.outfile);

/*
* Count up the number of flex label declarations and definitions,
* then allocate memory
* 6/95 te
*/
  for (label_flex->total = definition_count = 0; fgets (line, 100, flex_file);)
  {
    if (!strncmp (line, "name", 4))
      label_flex->total++;
    if (!strncmp (line, "definition", 10))
      definition_count++;
  }

  rewind (flex_file);

  definition_total = definition_count;

  ecalloc
  (
    (void **) &label_flex->member,
    label_flex->total + 1,
    sizeof (FLEX_MEMBER),
    "flexible bond labels",
    global.outfile
  );

  ecalloc
  (
    (void **) &label_flex->definition,
    definition_total,
    sizeof (NODE),
    "flexible bond label definitions",
    global.outfile
  );

  strcpy (label_flex->member[0].name, "null");

/*
* Read in the atom label definitions
* 6/95 te
*/
  for (label_flex->total = 1, definition_count = definition_total = 0;
    fgets (line, 100, flex_file);)
  {
/*
*   Process flexible label declaration
*   6/95 te
*/
    if (!strncmp (line, "name", 4))
    {
      if (label_flex->total > 1)
      {
        if (label_flex->member[label_flex->total - 1].drive_id == NEITHER)
          exit (fprintf (global.outfile,
            "ERROR get_flex_labels: Missing drive_id parameter in %s\n",
            label_flex->file_name));

        else if (label_flex->member[label_flex->total - 1].minimize_flag ==
          NEITHER)
          exit (fprintf (global.outfile,
            "ERROR get_flex_labels: Missing minimize parameter in %s\n",
            label_flex->file_name));

        else if (definition_total != 2)
          exit (fprintf (global.outfile,
            "ERROR get_flex_labels: Incorrect number of definitions in %s\n",
            label_flex->file_name));
      }

      label_flex->member[label_flex->total - 1].definition_total =
        definition_total;
      label_flex->member[label_flex->total].definition =
        &label_flex->definition[definition_count];
      definition_total = 0;

      if (sscanf
        (line, "%*s %s", label_flex->member[label_flex->total].name) < 1)
      {
        fprintf (global.outfile,
          "Incomplete label_flex->member declaration.\n");
        exit (EXIT_FAILURE);
      }

/*
*     Convert flex label name to lowercase
*     6/95 te
*/
      for (i = 0; i < strlen (label_flex->member[label_flex->total].name); i++)
        label_flex->member[label_flex->total].name[i] =
          (char) tolower (label_flex->member[label_flex->total].name[i]);

      label_flex->total++;

      label_flex->member[label_flex->total - 1].drive_id = NEITHER;
      label_flex->member[label_flex->total - 1].minimize_flag = NEITHER;
    }

/*
*   Process label_flex->member search identifier
*   10/95 te
*/
    else if (!strncmp (line, "drive_id", 6))
      sscanf
        (line, "%*s %d", &label_flex->member[label_flex->total - 1].drive_id);

/*
*   Process flex label minimize flag
*   10/95 te
*/
    else if (!strncmp (line, "minimize", 8))
      sscanf
        (line, "%*s %d",
        &label_flex->member[label_flex->total - 1].minimize_flag);

/*
*   Process flex label definition
*   6/95 te
*/
    else if (!strncmp (line, "definition", 10))
    {
      strtok (white_line (line), " ");

      if (!assign_node
        (&label_flex->definition[definition_count], TRUE))
      {
        fprintf (global.outfile, "Error assigning flex label definitions.\n");
        exit (EXIT_FAILURE);
      }

      definition_count++;
      definition_total++;
    }
  }

/*
* Update last flex label info also
* 6/95 te
*/
  label_flex->member[label_flex->total - 1].definition_total =
    definition_total;

  efclose (&flex_file);

/*
* Assign search torsion parameters if requested
* 10/96 te
*/
  if (label_flex->drive_flag)
    get_flex_search (label_flex);

/*
* Print out the flex label and their definitions
* 6/95 te

  if (global.output_volume == 'v')
  {
    fprintf (global.outfile, "\n____Flexible_Bond_Label_Definitions____\n\n");
    for (i = 0; i < label_flex->total; i++)
    {
      fprintf (global.outfile, "\n");
      fprintf (global.outfile, "%-20s%s\n", "name",
        label_flex->member[i].name);
      fprintf (global.outfile, "%-20s%d\n", "drive_id",
        label_flex->member[i].drive_id);
      fprintf (global.outfile, "%-20s%d\n", "minimize",
        label_flex->member[i].minimize_flag);

      fprintf (global.outfile, "%-20s%d\n", "positions",
        label_flex->member[i].torsion_total);

      fprintf (global.outfile, "%-20s", "torsions");
      for (j = 0; j < label_flex->member[i].torsion_total; j++)
        fprintf (global.outfile, "%g ", label_flex->member[i].torsion[j]);

      fprintf (global.outfile, "\n");
      for (j = 0; j < label_flex->member[i].definition_total; j++)
      {
        fprintf (global.outfile, "%-20s", "definition");
        print_node (&label_flex->member[i].definition[j], 0);
        fprintf (global.outfile, "\n");
      }

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

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

  return TRUE;
}


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

void free_flex_labels (LABEL_FLEX *label_flex)
{
  int i, j;

  for (i = 0; i < label_flex->total; i++)
  {
    efree ((void **) &label_flex->member[i].torsion);

    for (j = 0; j < label_flex->member[i].definition_total; j++)
      free_node (&label_flex->member[i].definition[j]);
  }

  efree ((void **) &label_flex->member);
  efree ((void **) &label_flex->definition);
}


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

Program to read in torsion.defn
11/95 Yax Sun
10/96 edited by Todd Ewing

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

int get_flex_search (LABEL_FLEX *label_flex)
{
  int i, j;			/* Counter variables */
  STRING200 line;		/* String used to compose output */
  char *token;
  FILE  *search_file;

  int   drive_id;
  int   torsion_total;
  float *torsion = NULL;	/* temporarily hold the torsion values */

  search_file = efopen (label_flex->search_file_name, "r", global.outfile);

  while (fgets (line, 200, search_file))
  {
    token = strtok (white_line (line), " ");

    if (!strcmp (token, "drive_id"))
    {
      if (token = strtok (NULL, " "))
        drive_id = atoi (token);

      else
        exit (fprintf (global.outfile,
          "ERROR get_flex_search: Search_id value not specified in %s\n",
          label_flex->search_file_name));

      if (!fgets (line, 200, search_file) ||
        !(token = strtok (white_line (line), " ")) ||
        strcmp (token, "positions"))
        exit (fprintf (global.outfile,
          "ERROR get_flex_search: Positions field doesn't follow Id in %s\n",
          label_flex->search_file_name));

      if (token = strtok (NULL, " "))
        torsion_total = atoi (token);

      else
        exit (fprintf (global.outfile,
          "ERROR get_flex_search: Postions value not specified in %s\n",
          label_flex->search_file_name));

      ecalloc
      (
        (void **) &torsion,
        torsion_total,
        sizeof (float),
        "torsion sample values",
        global.outfile
      );

      if (!fgets (line, 200, search_file) ||
        !(token = strtok (white_line (line), " ")) ||
        strcmp (token, "torsions"))
        exit (fprintf (global.outfile,
          "ERROR get_flex_search: Torsions doesn't follow Positions in %s\n",
          label_flex->search_file_name));

      for (i = 0; i < torsion_total; i++)
      {
        if (token = strtok (NULL, " "))
          torsion[i] = atof (token);

        else
          exit (fprintf (global.outfile,
            "ERROR get_flex_search: Insufficient number of torsions in %s\n",
            label_flex->search_file_name));
      }

      for (i = 0; i < label_flex->total; i++)
      {
        if (label_flex->member[i].drive_id == drive_id)
        {
          label_flex->member[i].torsion_total = torsion_total;

          ecalloc
          (
            (void **) &label_flex->member[i].torsion,
            label_flex->member[i].torsion_total,
            sizeof (float),
            "stored torsion sample values",
            global.outfile
          );

          for (j = 0; j < label_flex->member[i].torsion_total; j++)
            label_flex->member[i].torsion[j] = torsion[j];

/*
          fprintf (global.outfile,
            "search %d, %s\n", drive_id, label_flex->member[i].name);
*/
        }
      }

      torsion_total = 0;
      efree ((void **) &torsion);
    }
  }

  for (i = 1; i < label_flex->total; i++)
    if (label_flex->member[i].torsion_total < 1)
      exit (fprintf (global.outfile,
        "ERROR get_flex_search: Missing torsion parameters in %s\n",
        label_flex->search_file_name));

  fclose (search_file);
  return TRUE;
}


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

void assign_flex_labels
(
  LABEL_FLEX	*label_flex,
  MOLECULE	*molecule
)
{
  int i, j;
  int bond_id, bond_found;
  int flex_id;
  float conf_total;
  static MOLECULE temporary = {0};

  void detect_rings (MOLECULE *);

  if (!label_flex->init_flag)
    get_flex_labels (label_flex);

  reset_molecule (&temporary);

/*
* Make a copy of the original set of torsion bonds
* 10/96 te
*/
  copy_torsions (&temporary, molecule);

/*
* Flag all ring bonds
* 10/96 te
*/
  detect_rings (molecule);

/*
* Identify all bonds with flex labels
* 10/96 te
*/
  for (i = molecule->total.torsions = 0; i < molecule->total.bonds; i++)
  {
    molecule->bond[i].flex_id = 0;

    if (molecule->bond[i].ring_flag == TRUE)
      continue;

    if ((molecule->atom[molecule->bond[i].origin].neighbor_total <= 1) ||
      (molecule->atom[molecule->bond[i].target].neighbor_total <= 1))
      continue;

    for (j = 1; j < label_flex->total; j++)
    {
      if
      (
        check_atom
          (molecule, molecule->bond[i].origin,
          &label_flex->member[j].definition[0]) &&
        check_atom
          (molecule, molecule->bond[i].target,
          &label_flex->member[j].definition[1])
      )
        molecule->bond[i].flex_id = j;

      else if
      (
        check_atom
          (molecule, molecule->bond[i].origin,
          &label_flex->member[j].definition[1]) &&
        check_atom
          (molecule, molecule->bond[i].target,
          &label_flex->member[j].definition[0])
      )
        molecule->bond[i].flex_id = j;
    }

    if (molecule->bond[i].flex_id)
      molecule->total.torsions++;
  }

/*
* Check for sufficient space to store new and old torsions
* 10/96 te
*/
  if ((molecule->total.torsions + temporary.total.torsions) >
    molecule->max.torsions)
  {
    molecule->total.torsions += temporary.total.torsions;
    reallocate_torsions (molecule);
  }

  else
    reset_torsions (molecule);

/*
* Store flexible torsions
* 10/96 te
*/
  for (i = molecule->total.torsions = 0; i < molecule->total.bonds; i++)
    if (molecule->bond[i].flex_id)
    {
      molecule->torsion[molecule->total.torsions].bond_id = i;
      molecule->torsion[molecule->total.torsions].flex_id =
        molecule->bond[i].flex_id;
      molecule->total.torsions++;
    }

/*
* Append original torsions
* 10/96 te
*/
  for (i = 0; i < temporary.total.torsions; i++)
  {
    for (j = 0, bond_found = FALSE; j < molecule->total.torsions; j++)
      if (temporary.torsion[i].bond_id == molecule->torsion[j].bond_id)
      {
        molecule->torsion[j].target_angle = temporary.torsion[i].target_angle;
        bond_found = TRUE;
      }

    if (!bond_found)
      copy_torsion
        (&molecule->torsion[molecule->total.torsions++],
        &temporary.torsion[i]);
  }

/*
* Compute current torsion angles and conformation total
* 10/96 te
*/
  get_torsion_neighbors (molecule);

  for (i = 0, conf_total = 1.0; i < molecule->total.torsions; i++)
  {
    molecule->torsion[i].current_angle =
      molecule->torsion[i].target_angle =
      compute_torsion (molecule, i);

    molecule->torsion[i].periph_flag =
      check_peripheral_torsion (molecule, i);

    bond_id = molecule->torsion[i].bond_id;

    if (label_flex->drive_flag)
    {
      flex_id = molecule->bond[bond_id].flex_id;

      if (flex_id > 0)
      {
        if (label_flex->member[flex_id].torsion_total < 1)
          exit (fprintf (global.outfile,
            "ERROR assign_flex_labels: bond with no positions detected\n"));

        conf_total *= (float) label_flex->member[flex_id].torsion_total;
      }
    }
  }

/*
* Print out flex label assignments
* 10/95 te
*/
  if (global.output_volume == 'v')
  {
    fprintf (global.outfile, "____Flexible_Bond_Label_Assignments____\n");
    fprintf (global.outfile, "%s:\n\n", molecule->info.name);
    fprintf (global.outfile,
      "%-4s %-4s | %-4s %-4s %-4s %-4s | %7s %3s %3s %s\n",
      "tors", "bond", "orig", "orig", "targ", "targ", "input", "pos", "drv", "flex");
    fprintf (global.outfile,
      "%-4s %-4s | %-4s %-4s %-4s %-4s | %7s %3s %3s %s\n",
      " id", " id", "nhbr", "", "", "nhbr", "angle", "num", "id", "label");
    fprintf (global.outfile,
      "-----------------------------------------"
      "---------------------------\n");

    for (i = 0; i < molecule->total.torsions; i++)
    {
      bond_id = molecule->torsion[i].bond_id;

      fprintf
      (
        global.outfile,
        "%-4d %-4d | %-4s %-4s %-4s %-4s | %7.2f %3d %3d %s\n",
        i + 1, bond_id + 1,
        molecule->atom[molecule->torsion[i].origin_neighbor].name,
        molecule->atom[molecule->torsion[i].origin].name,
        molecule->atom[molecule->torsion[i].target].name,
        molecule->atom[molecule->torsion[i].target_neighbor].name,
        molecule->torsion[i].current_angle / PI * 180,
        label_flex->member[molecule->bond[bond_id].flex_id].torsion_total,
        molecule->bond[bond_id].flex_id ?
          label_flex->member[molecule->bond[bond_id].flex_id].drive_id :
          0,
        molecule->bond[bond_id].flex_id ?
          label_flex->member[molecule->bond[bond_id].flex_id].name :