io_pdb.c

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

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

/*
Written by Todd Ewing
10/95
*/
#include "define.h"
#include "utility.h"
#include "mol.h"
#include "global.h"
#include "io_pdb.h"

typedef struct link_atom
{
  ATOM atom;
  XYZ coord;
  struct link_atom *previous, *next;
} LINK_ATOM;

typedef struct link_bond
{
  BOND bond;
  struct link_bond *previous, *next;
} LINK_BOND;

typedef struct link_subst
{
  SUBST subst;
  struct link_subst *previous, *next;
} LINK_SUBST;


int read_pdb
(
  MOLECULE *molecule,
  FILE_NAME molecule_file_name,
  FILE *molecule_file
)
{
  int i, j, origin, target;
  char buff[200], buff2[200];
  STRING20 temp, subst_name;
  int subst_number;
  int found_atom;
  int new_substructure;
  LINK_ATOM *first_atom, *previous_atom, *current_atom;
  LINK_BOND *first_bond, *previous_bond, *current_bond;
  LINK_SUBST *first_subst, *previous_subst, *current_subst;

  void assign_atom_type (ATOM *);

/*
* Initialize molecule info
* 6/95 te
*/
  vstrcpy (&molecule->info.name, "****");
  vstrcpy (&molecule->info.status_bits, "****");
  vstrcpy (&molecule->info.comment, "****");

/*
* Read in molecule info if found before ATOM records
* 6/95 te
*/
  if (!fgets (buff, 199, molecule_file))
    return EOF;

  while (strncmp (buff, "ATOM", 4) && strncmp (buff, "HETATM", 6))
  {
    if (!strncmp (buff, "HEADER", 6))
    {
      sscanf (buff, "%*s %s", buff2);
      vstrcpy (&molecule->info.name, buff2);
    }

    if (!strncmp (buff, "COMPND", 6))
    {
      for (i = 6; (i < strlen (buff)) && isspace (buff[i]); i++);
      vstrcpy (&molecule->info.comment, strip_newline (&buff[i]));
      if (strchr (molecule->info.comment, ' '))
        strcpy (strchr (molecule->info.comment, ' '), "");
    }

    if (!fgets (buff, 199, molecule_file))
      return EOF;
  }

/*
* Initialize linked lists and counters
* 6/95 te
*/
  first_atom = current_atom = previous_atom = NULL;
  first_subst = current_subst = previous_subst = NULL;
  first_bond = current_bond = previous_bond = NULL;

/*
* Read in atom information and store in a linked list
* 6/95 te
*/
  while (!strncmp (buff, "ATOM", 4) || !strncmp (buff, "HETATM", 6))
  {
/*
*   Allocate memory for this ATOM link
*   6/95 te
*/
    previous_atom = current_atom;
    current_atom = NULL;

    ecalloc
    (
      (void **) &current_atom,
      1,
      sizeof (LINK_ATOM),
      "linked atom list",
      global.outfile
    );

/*
*   Update pointers between ATOM links, initialize first link
*   6/95 te
*/
    current_atom->previous = previous_atom;

    if (previous_atom)
      previous_atom->next = current_atom;

    else
      first_atom = current_atom;

/*
*   Deposit atom info in this ATOM link
*   6/95 te
*/
    memset (temp, 0, sizeof (STRING20));
    current_atom->atom.number = atoi (strncpy (temp, &buff[6], 5));

    memset (buff2, 0, 199);
    strncpy (buff2, &buff[12], 4);
    sscanf (buff2, "%s", temp);
    vstrcpy (&current_atom->atom.name, temp);

    memset (temp, 0, sizeof (STRING20));
    strncpy (subst_name, &buff[17], 3);
    subst_name[3] = 0;

    memset (temp, 0, sizeof (STRING20));
    subst_number = atoi (strncpy (temp, &buff[22], 4));

    memset (temp, 0, sizeof (STRING20));
    current_atom->coord[0] = atof (strncpy (temp, &buff[30], 8));

    memset (temp, 0, sizeof (STRING20));
    current_atom->coord[1] = atof (strncpy (temp, &buff[38], 8));

    memset (temp, 0, sizeof (STRING20));
    current_atom->coord[2] = atof (strncpy (temp, &buff[46], 8));

    if (!strcmp (strrchr (molecule_file_name, '.'), ".xpdb"))
    {
      memset (temp, 0, sizeof (STRING20));
      current_atom->atom.charge = atof (strncpy (temp, &buff[55], 8));

      memset (temp, 0, sizeof (STRING20));
      strncpy (current_atom->atom.type, &buff[71], 5);
    }

    else
    {
      current_atom->atom.charge = 0.0;
      assign_atom_type (&current_atom->atom);
    }

/*
*   Store substructure information if the current residue is new
*   6/95 te
*/
    if (!previous_atom || !current_subst ||
      (subst_number != current_subst->subst.number))
    {
/*
*     Check to see if this substructure has been seen before
*     6/95 te
*/
      for (i = 0, current_subst = first_subst, new_substructure = TRUE;
        current_subst != NULL;
        i++, previous_subst = current_subst,
        current_subst = current_subst->next)
      {
        if (current_subst->subst.number == subst_number)
        {
          new_substructure = FALSE;
          current_atom->atom.subst_id = i;
        }
      }

/*
*     If the current atom is a member of a new residue, then add a new link
*     10/95 te
*/
      if (new_substructure)
      {
        current_subst = NULL;

        ecalloc
        (
          (void **) &current_subst,
          1,
          sizeof (LINK_SUBST),
          "linked subst list",
          global.outfile
        );

/*
*       Update pointers between SUBST links, initialize first link
*       6/95 te
*/
        current_subst->previous = previous_subst;

        if (previous_subst)
          previous_subst->next = current_subst;
        else
          first_subst = current_subst;

/*
*       Deposit substructure info in this SUBST link
*       6/95 te
*/
        current_subst->subst.number = subst_number;
        vstrcpy (&current_subst->subst.name, subst_name);
        current_subst->subst.root_atom = molecule->total.atoms;

        current_atom->atom.subst_id = molecule->total.substs++;
      }

/*
*     Otherwise, set current subst pointer to the last link added
*     10/95 te
*/
      else
        current_subst = previous_subst;
    }

    else
      current_atom->atom.subst_id = molecule->total.substs - 1;

/*
*   Update atom info
*   10/95 te
*/
    molecule->total.atoms++;

/*
*   Read in next line of file
*   6/95 te
*/
    if (!fgets (buff, 199, molecule_file))
      break;
  }

/*
* Advance to TER record, but stop if CONECT records are found
* 6/95 te
*/
  while (strncmp (buff, "TER", 3))
  {
    if (!strncmp (buff, "CONECT", 6))
      break;

    if (!fgets (buff, 199, molecule_file))
      break;
  }

/*
* Read in connectivity information
* 6/95 te
*/
  while (!strncmp (buff, "CONECT", 6))
  {
    memset (temp, 0, sizeof (STRING20));
    origin = atoi (strncpy (temp, &buff[6], 5));

/*
*   Scan the data in the CONECT record
*   6/95 te
*/
    for (i = 1; target = atoi (strncpy (temp, &buff[i * 5 + 6], 5)); i++)
    {
      if (target > origin)
      {
/*
*       Allocate memory for this BOND link
*       6/95 te
*/
        previous_bond = current_bond;
        current_bond = NULL;

        ecalloc
        (
          (void **) &current_bond,
          1,
          sizeof (LINK_BOND),
          "linked bond list",
          global.outfile
        );

/*
*       Update pointers between BOND links, initialize first link
*       6/95 te
*/
        current_bond->previous = previous_bond;

        if (previous_bond)
          previous_bond->next = current_bond;

        else
          first_bond = current_bond;

/*