io_mol2.c

来自「最经典的分子对结软件」· C语言 代码 · 共 992 行 · 第 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_mol2.h"

int read_mol2
(
  MOLECULE *molecule,
  FILE_NAME molecule_file_name,
  FILE *molecule_file
)
{
  int i, j;			/* Iteration variables */
  int token_error = FALSE;	/* Flag for error reading tokens */
  char *line = NULL;		/* Line of data */
  char *token;			/* Pointer to identify each token */
  int return_value = TRUE;	/* Value to return upon completion */

/*
* Find the next molecule data record
* 6/96 te
*/
  if (find_record (&line, "@<TRIPOS>MOLECULE", molecule_file) == NULL)
  {
    return_value = EOF;
    goto terminate;
  }

/*
* Record molecule name.  If no name, then record default name.
* 10/96 te
*/
  if (vfgets (&line, molecule_file) == NULL)
  {
    return_value = EOF;
    goto terminate;
  }

  for (token = strtok (white_line (line), " "); token;
    token = strtok (NULL, " "))
  {
    if (molecule->info.name)
      vstrcat (&molecule->info.name, " ");
    vstrcat (&molecule->info.name, token);
  }
    
  if (!molecule->info.name || !strcmp (molecule->info.name, ""))
    vstrcpy (&molecule->info.name, "****");

/*
* Record the total molecule components and types
* 10/96 te
*/
  if (vfgets (&line, molecule_file) == NULL)
  {
    return_value = EOF;
    goto terminate;
  }

  if (token = strtok (white_line (line), " "))
    molecule->total.atoms = atoi (token);

  else
  {
    fprintf (global.outfile,
      "WARNING read_mol2: incorrect MOLECULE record for %s in %s\n",
      molecule->info.name, molecule_file_name);
    return_value = NULL;
    goto terminate;
  }

  if (token = strtok (NULL, " "))
  {
    molecule->total.bonds = atoi (token);

    if (token = strtok (NULL, " "))
    {
      molecule->total.substs = atoi (token);

      token = strtok (NULL, " ");
        /* Skip number features data record */

      if (token = strtok (NULL, " "))
        molecule->total.sets = atoi (token);

      else
        molecule->total.sets = 0;
    }

    else
    {
      molecule->total.substs = 0;
      molecule->total.sets = 0;
    }
  }

  else
  {
    molecule->total.bonds = 0;
    molecule->total.substs = 0;
    molecule->total.sets = 0;
  }

  if (vfgets (&line, molecule_file) == NULL)
  {
    return_value = EOF;
    goto terminate;
  }

  vstrcpy (&molecule->info.molecule_type, strtok (white_line (line), " "));

  if (vfgets (&line, molecule_file) == NULL)
  {
    return_value = EOF;
    goto terminate;
  }

  vstrcpy (&molecule->info.charge_type, strtok (white_line (line), " "));

  if (vfgets (&line, molecule_file) == NULL)
  {
    return_value = EOF;
    goto terminate;
  }

  if (strncmp (line, "@<TRIPOS>", 9))
  {
    vstrcpy (&molecule->info.status_bits, strtok (white_line (line), " "));

/*
*   Record molecule comment.  If no comment was read,
*   then record default comment.
*   10/96 te
*/
    if (vfgets (&line, molecule_file) == NULL)
  {
    return_value = EOF;
    goto terminate;
  }

    for (token = strtok (white_line (line), " "); token;
      token = strtok (NULL, " "))
    {
      if (molecule->info.comment)
        vstrcat (&molecule->info.comment, " ");
      vstrcat (&molecule->info.comment, token);
    }
    
    if (!molecule->info.comment || !strcmp (molecule->info.comment, ""))
      vstrcpy (&molecule->info.comment, "****");
  }

  else
  {
    vstrcpy (&molecule->info.status_bits, "****");
    vstrcpy (&molecule->info.comment, "****");
    fseek (molecule_file, (long) -strlen (line), SEEK_CUR);
  }

/*
* Allocate space for molecule components
* 2/96 te
*/
  reallocate_atoms (molecule);
  reallocate_bonds (molecule);
  reallocate_substs (molecule);
  reallocate_sets (molecule);

/*
* Read in atoms
* 6/95 te
*/
  if (find_record (&line, "@<TRIPOS>ATOM", molecule_file) == NULL)
  {
    fprintf (global.outfile,
      "WARNING read_mol2: Unable to find ATOM record for %s in %s\n",
      molecule->info.name, molecule_file_name);
    return_value = NULL;
    goto terminate;
  }

  for (i = 0; i < molecule->total.atoms; i++)
  {
    if (vfgets (&line, molecule_file) == NULL)
    {
      fprintf (global.outfile,
        "WARNING read_mol2: Incomplete ATOM record for %s in %s\n",
        molecule->info.name, molecule_file_name);
      return_value = NULL;
      goto terminate;
    }

/*
*   Parse each data line into space-separated tokens and check that each
*   token is found.
*   9/95 te
*/
    molecule->atom[i].number = i + 1;

    if (token = strtok (white_line (line), " "));
      /* Skip the atom id field */
    else
      token_error = TRUE;

    if (!token_error && (token = strtok (NULL, " ")))
      vstrcpy (&molecule->atom[i].name, token);
    else
      token_error = TRUE;

    if (!token_error && (token = strtok (NULL, " ")))
      molecule->coord[i][0] = atof (token);
    else
      token_error = TRUE;

    if (!token_error && (token = strtok (NULL, " ")))
      molecule->coord[i][1] = atof (token);
    else
      token_error = TRUE;

    if (!token_error && (token = strtok (NULL, " ")))
      molecule->coord[i][2] = atof (token);
    else
      token_error = TRUE;

    if (!token_error && (token = strtok (NULL, " ")))
      vstrcpy (&molecule->atom[i].type, token);
    else
      token_error = TRUE;

    if (!token_error && (token = strtok (NULL, " ")))
      molecule->atom[i].subst_id = atoi (token) - 1;
    else
      token_error = TRUE;

    if ((molecule->atom[i].subst_id < 0) ||
      (molecule->atom[i].subst_id >= molecule->total.substs))
    {
      if (molecule->total.substs == 1)
        molecule->atom[i].subst_id = 0;

      else
      {
        fprintf (global.outfile,
          "WARNING read_mol2: "
          "Improper ATOM record substructure id for %s in %s.\n",
          molecule->info.name, molecule_file_name);
        return_value = NULL;
        goto terminate;
      }
    }

    if (!token_error && (token = strtok (NULL, " ")));
      /* Skip the substructure name field */
    else
      token_error = TRUE;

    if (!token_error && (token = strtok (NULL, " ")))
      molecule->atom[i].charge = atof (token);
    else
      token_error = TRUE;

    if (token_error)
    {
      fprintf (global.outfile,
        "WARNING read_mol2: Error reading ATOM record for %s in %s.\n",
        molecule->info.name, molecule_file_name);
      return_value = NULL;
      goto terminate;
    }
  }

/*
* Read in bonds
* 6/95 te
*/
  if (molecule->total.bonds > 0)
  {
    if (find_record (&line, "@<TRIPOS>BOND", molecule_file) == NULL)
    {
      fprintf (global.outfile,
        "WARNING read_mol2: Unable to fine BOND record for %s in %s\n",
        molecule->info.name, molecule_file_name);
      return_value = NULL;
      goto terminate;
    }

    for (i = 0; i < molecule->total.bonds; i++)
    {
      if (vfgets (&line, molecule_file) == NULL)
      {
        fprintf (global.outfile,
          "WARNING read_mol2: Incomplete BOND record for %s in %s\n",
          molecule->info.name, molecule_file_name);
        return_value = NULL;
        goto terminate;
      }

/*
*     Parse each data line into space-separated tokens and check that each
*     token is found.
*     9/95 te
*/
      if (token = strtok (white_line (line), " "));
        /* Skip the bond id field */
      else
        token_error = TRUE;

      if (!token_error && (token = strtok (NULL, " ")))
        molecule->bond[i].origin = atoi (token) - 1;
      else
        token_error = TRUE;

      if (!token_error && (token = strtok (NULL, " ")))
        molecule->bond[i].target = atoi (token) - 1;
      else
        token_error = TRUE;

      if (!token_error && (token = strtok (NULL, " ")))
        vstrcpy (&molecule->bond[i].type, strip_char (token, '\n'));
      else
        token_error = TRUE;

      if (token_error)
      {
        fprintf (global.outfile,
          "WARNING read_mol2: Error in BOND record of %s in %s\n",
          molecule->info.name, molecule_file_name);
        return_value = NULL;
        goto terminate;
      }

      if ((molecule->bond[i].origin >= molecule->total.atoms) ||
        (molecule->bond[i].origin < 0) ||
        (molecule->bond[i].target >= molecule->total.atoms) ||
        (molecule->bond[i].target < 0))
      {
        fprintf (global.outfile,
          "WARNING read_mol2: Improper BOND origin/target for %s in %s.\n",
          molecule->info.name, molecule_file_name);
        return_value = NULL;
        goto terminate;
      }
    }
  }

/*
* Read in substructure information
* 9/95 te
*/
  if (molecule->total.substs > 0)
  {
    if (find_record (&line, "@<TRIPOS>SUBSTRUCTURE", molecule_file) == NULL)
    {
      fprintf (global.outfile,
        "WARNING read_mol2: No SUBSTRUCTURE record for %s in %s\n",
        molecule->info.name, molecule_file_name);
      return_value = NULL;
      goto terminate;
    }

    for (i = 0; i < molecule->total.substs; i++)
    {
      if (vfgets (&line, molecule_file) == NULL)
      {
        fprintf (global.outfile,
          "WARNING read_mol2: Incomplete SUBSTRUCTURE record for %s in %s\n",
          molecule->info.name, molecule_file_name);
        return_value = NULL;
        goto terminate;
      }
/*
*     Parse each data line into space-separated tokens and check that each
*     token is found.
*     9/95 te
*/
      molecule->subst[i].number = i + 1;

      if (token = strtok (white_line (line), " "));
        /* Skip the substructure id field */
      else
        token_error = TRUE;

      if (!token_error && (token = strtok (NULL, " ")))
        vstrcpy (&molecule->subst[i].name, token);
      else
        token_error = TRUE;

      if (!token_error && (token = strtok (NULL, " ")))
        molecule->subst[i].root_atom = atoi (token) - 1;
      else
        token_error = TRUE;

      if (!token_error && (token = strtok (NULL, " ")))
      {
        vstrcpy (&molecule->subst[i].type, strip_char (token, '\n'));

        if (token = strtok (NULL, " "))
        {
          molecule->subst[i].dict_type = atoi (token) - 1;

          if (token = strtok (NULL, " "))
          {
            vstrcpy (&molecule->subst[i].chain, strip_char (token, '\n'));

            if (token = strtok (NULL, " "))
            {
              vstrcpy (&molecule->subst[i].sub_type, strip_char (token, '\n'));

              if (token = strtok (NULL, " "))
              {
                molecule->subst[i].inter_bonds = atoi (token);

                if (token = strtok (NULL, " "))
                  vstrcpy (&molecule->subst[i].status, strip_char (token, '\n'));
              }
            }
          }
        }
      }

      if (token_error)
      {
        fprintf (global.outfile,
          "WARNING read_mol2: Error in SUBSTRUCTURE record for %s in %s\n",
          molecule->info.name, molecule_file_name);
        return_value = NULL;
        goto terminate;
      }
    }

/*
*   Update substructure information if more that one substructure is present.
*   (Assume molecule is a macromolecule.)  Extract out the residue number
*   from the substructure name only if the beginning of the name is the
*   same as the sub_type.
*   10/95 te
*/
    if (molecule->total.substs > 1)
    {
      for (i = 0, token_error = TRUE; i < molecule->total.substs; i++)
        if ((strlen (molecule->subst[i].name) <= 3) ||
          !isdigit (molecule->subst[i].name[3]) ||
          (atoi (&molecule->subst[i].name[3]) <= 0))
          token_error = FALSE;
     
      if (token_error)
        for (i = 0; i < molecule->total.substs; i++)
          molecule->subst[i].number = atoi
            (&molecule->subst[i].name[3]);
    }
  }

  else
  {
    molecule->total.substs = 1;
    reallocate_substs (molecule);

    molecule->subst[0].number = 1;
    vstrcpy (&molecule->subst[0].name, "****");
    molecule->subst[0].root_atom = 0;

    for (i = 0; i < molecule->total.atoms; i++)
      molecule->atom[i].subst_id = 0;
  }

/*
* Read in set information
* 10/96 te
*/
  if (molecule->total.sets > 0)
  {
    token_error = FALSE;

    if (find_record (&line, "@<TRIPOS>SET", molecule_file) == NULL)
    {
      fprintf (global.outfile,
        "WARNING read_mol2: No SET record for %s in %s\n",
        molecule->info.name, molecule_file_name);
      return_value = NULL;
      goto terminate;
    }

    for (i = 0; i < molecule->total.sets; i++)
    {
      if (vfgets (&line, molecule_file) == NULL)
      {
        fprintf (global.outfile,