dock.c

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

/***************************************************************************\
*                                                                           *
*                                                                           *
*               UUUUUUUUU     CCCCCCC     SSSSSSS    FF/   FFF/             *
*               UU/    UU/  CC/    CC/  SS/    SS/  FF/ FFF/                *
*              UU/    UU/  CC/    CC/  SS/         FFFFF/                   *
*             UU/    UU/  CC/    CC/  SS/         FF/ FF\                   *
*            UU/    UU/  CC/    CC/  SS/    SS/  FF/   FF\                  *
*          UUUUUUUUU/    CCCCCCC/    SSSSSSS/   FF/     FF\                 *
*                                                                           *
*                                                                           *
*       Copyright (C) 1997 Regents of the University of California          *
*                         All Rights Reserved.                              *
*                                                                           *
*     This program implements the docking algorithm of I.D. Kuntz,          *
*     J Mol Biol 161, 269-288, 1982.                                        *
*                                                                           *
*     Versions 1.0 and 1.1 of the dock program were based extensively on    *
*     the work of Robert Sheridan, Renee DesJarlais, and Tack Kuntz.        *
*                                                                           *
*     Version 2.0, macrodock, is based on the work of Brian Shoichet        *
*     and Tack Kuntz, with help from Dale Bodian.                           *
*                                                                           *
*     Version 3.0, chemdock, is based on the work of Elaine Meng, Brian     *
*     Shoichet, and Tack Kuntz.                                             *
*                                                                           *
*     Version 3.5 was produced from version 3.0 by Mike Connolly and        *
*     Dan Gschwend.                                                         *
*                                                                           *
*     Version 4.0 is based on the work of Todd Ewing and Tack Kuntz.        *
*                                                                           *
*                                                                           *
\***************************************************************************/

#include <time.h>
#include "define.h"
#include "global.h"
#include "utility.h"
#include "mol.h"
#include "search.h"
#include "dock.h"
#include "label.h"
#include "screen.h"
#include "score.h"
#include "score_dock.h"
#include "flex.h"
#include "match.h"
#include "orient.h"
#include "io.h"
#include "io_grid.h"
#include "io_ligand.h"
#include "mol_prep.h"
#include "parm_dock.h"
#include "rank.h"

GLOBAL global = {0};		/* See global.h for global variables */

main (int argc, char *argv[])
{
/*
* General variables (see header files for globally defined variables
* 6/95 te
*/
  int i;			/* Counter variables */
  float time = 0.0;		/* Total amount of elapsed cpu time */
  int write_flag;		/* Flag to write molecule when scored */

/*
* Variables for multiple ligands mode
* 6/95 te
*/
  int lig_flag;			/* Flag for whether a ligand was read */
  int lig_total = 0;		/* Number of ligands read */
  int lig_proc = 0;		/* Number of ligands processed */
  int lig_skip = 0;		/* Number of ligands skipped */

/*
* File pointers used when checking for the presence of control files
* 6/95 te
*/
  FILE *dump = NULL;
  FILE *quit = NULL;

/*
* Data structures containing docking parameters
* 2/96 te
*/
  DOCK dock = {0};		/* Docking data structure */
  ORIENT orient = {0};		/* Orientation data structure */
  SCORE score = {0};		/* Scoring data structure */
  LABEL label = {0};		/* Labeling data structure */

  LIST best_anchors = {0};	/* List of best anchor orientations */
  LIST best_orients = {0};	/* List of best orientations */
  LIST best_ligands = {0};	/* List of best molecules */

  MOLECULE mol_ref = {0};	/* Reference ligand data (from file) */
  MOLECULE mol_init = {0};	/* Initial ligand data */
  MOLECULE mol_conf = {0};	/* Ligand conformation data */
  MOLECULE mol_ori = {0};	/* Oriented ligand data */
  MOLECULE mol_score = {0};	/* Minimized ligand data */
  MOLECULE mol_out = {0};	/* Output ligand data */

/*
* Functions used in the main routine that also reside in this file.
* Other functions are declared in the header files.
* 6/95 te
*/
  void write_program_header (void);
  void set_memory_limit (void);
  float elapsed_time (float *);
  void initialize_performance (DOCK *, SCORE *);
  void report_performance (DOCK *, SCORE *, int);

/*
* Process command line arguments
* 6/95 te
*/
  process_commands (&dock, argc, argv);

/*
* Output program header, set memory ceiling, and initialize clock
* 6/95 te
*/
  write_program_header ();
  set_memory_limit ();

  if (dock.performance_flag)
    initialize_performance (&dock, &score);

  dock.total_time = -elapsed_time (NULL);

/*
* Read in dock parameters
* 6/95 te
*/
  get_parameters (&dock, &orient, &score, &label);

/*
* Initialize best molecule lists
* 12/96 te
*/
  allocate_lists (&score, &best_anchors, dock.rank_anchor_total, FALSE);
  allocate_lists (&score, &best_orients, dock.rank_orient_total, FALSE);
  allocate_lists (&score, &best_ligands, dock.rank_ligand_total, TRUE);

/*
* Open ligand input files 
* 6/95 te
*/
  if (!dock.parallel.flag || dock.parallel.server)
    dock.ligand_file = efopen (dock.ligand_file_name, "r", global.outfile);

/*
* Open ligand output files 
* 6/95 te
*/
  if (!dock.rank_ligands && !dock.parallel.server)
    for (i = 0; i < SCORE_TOTAL; i++)
      if (score.type[i].flag)
        score.type[i].file =
          efopen (score.type[i].file_name, "w", global.outfile);

  write_flag =
    (dock.write_orients && !dock.rank_orients) ||
    (!dock.multiple_orients && dock.multiple_ligands && !dock.rank_ligands);

/*
* If a restart has been signalled, then read in previous data
* 11/96 te
*/
  if (dock.restart)
  {
    dock.restart = FALSE;
    fprintf (global.outfile,
      "Reading restart information from disk.\n");

    if (read_restartinfo
    (
      &dock,
      &score,
      &best_ligands,
      &lig_total,
      &lig_proc,
      &lig_skip,
      &time
    ))
    {
/*
*     Reset time to what it was in previous run
*     11/95 te
*/
      elapsed_time (&time);
    }

    else
      fprintf (global.outfile,
        "WARNING main: Restart information not found."
        " Normal docking to be performed.\n");

    dock.total_time =
      dock.read_time =
      -elapsed_time (NULL);
  }

/*
* Otherwise, skip initial set of molecules
* 11/96 te
*/
  else
  {
    dock.read_time -= elapsed_time (NULL);

    for
    (
      i = 0;
      (i < dock.initial_skip) &&
        (get_ligand
        (
          &dock, &score, &label,
          &mol_ref, &mol_init,
          FALSE, FALSE, FALSE
        ) != EOF);
      lig_total = lig_skip = ++i
    );
  }

/*
* Loop over all ligands
* 6/95 te
*/
  for (;
    (!lig_total ||
      (dock.multiple_ligands && (lig_total < dock.max_ligands))) &&
    ((lig_flag = get_ligand
      (
        &dock, &score, &label,
        &mol_ref, &mol_init,
        score.flag || label.flex.flag,
        label.chemical.flag,
        label.vdw.flag
      )) != EOF);
    lig_total++)
  {
/*
*   Check to see if any ligands were read from input file
*   8/95 te
*/
    if (lig_flag == TRUE) 
    {
      if (global.output_volume == 'v')
      {
        fprintf (global.outfile, "Processing %s\n", mol_ref.info.name);
        fflush (global.outfile);
      }

/*
*     Perform chemical screening functions
*     11/96 te
*/
      if (label.chemical.screen.process_flag)
      {
        dock.screen_time -= elapsed_time (NULL);

        if (check_screen (&orient.match, &label, &mol_ref, lig_proc))
        {
          if (dock.rank_ligands)
            update_list (&best_ligands, &mol_ref);

          else
            write_ligand
            (
              &dock,
              &score,
              &mol_ref,
              score.type[NONE].file_name,
              score.type[NONE].file
            );

          fprintf (global.outfile, "+");
        }

        else
          fprintf (global.outfile, ".");

        if (++lig_proc % 50 == 0)
          fprintf (global.outfile, "| %d mols, %.4g secs\n",
            lig_proc, elapsed_time (NULL));

        fflush (global.outfile);
        dock.screen_time += elapsed_time (NULL);
        continue;
      }

/*
*     Initialize ligand data
*     10/96 te
*/
      if (dock.multiple_ligands && (score.flag || orient.flag))
        mol_ref.info.input_id =
          mol_init.info.input_id = lig_total + 1;

      dock.orient_total = dock.score_total = dock.conform_total = 0;
      reset_lists (&score, &best_orients);

/*
*     Loop over all anchors
*     12/96 te
*/
      dock.conform_time -= elapsed_time (NULL);

      for
      (
        dock.anchor_total = 0;
        get_anchor
        (
          &label.flex,
          &mol_init,
          &mol_conf,
          dock.anchor_total
        ) != EOF;
        dock.anchor_total++
      )
      {
        reset_lists (&score, &best_anchors);

/*
*       Loop over all anchor conformations
*       12/96 te
*/
        for
        (
          dock.conform_subtotal = 0;
          get_anchor_conformation
          (
            &label,
            &score,
            &mol_init,
            &mol_conf,
            dock.conform_subtotal
          ) != EOF;
          dock.conform_subtotal++, dock.conform_total++
        )
        {
/*
*         Update chemical screen keys in ligand
*         11/96 te
*/
          if (label.chemical.screen.construct_flag)
          {
            update_keys (&label.chemical, &mol_conf, dock.conform_subtotal);
            continue;
          }

          else if (score.flag || orient.flag)
          {
/*
*           Initialize ligand data structures
*           10/96 te
*/
            if (dock.conform_subtotal == 0)
            {
              copy_molecule (&mol_ori, &mol_conf);
              copy_molecule (&mol_score, &mol_conf);
            }

/*
*           Update ligand orientation data structure
*           10/96 te
*/
            else
            {
              copy_transform (&mol_ori, &mol_conf);
              copy_score (&mol_ori.score, &mol_conf.score);
              copy_coords (&mol_ori, &mol_conf);
              copy_torsions (&mol_ori, &mol_conf);
              copy_segments (&mol_ori, &mol_conf);
              copy_layers (&mol_ori, &mol_conf);
            }

/*
*           Loop over all orientations
*           11/96 te
*/
            dock.orient_time -= elapsed_time (NULL);

            for
            (
              dock.orient_subtotal = dock.score_subtotal = 0;
              (get_orientation
              (
                &dock, &orient, &label,
                &mol_ref, &mol_conf, &mol_ori,
                lig_proc,
                dock.conform_subtotal,
                dock.orient_subtotal) != EOF) &&
              ((score.flag && (dock.score_subtotal < orient.max) &&
                (dock.orient_subtotal < 100 * orient.max)) ||
              (!score.flag && (dock.orient_subtotal < orient.max)));
              dock.orient_subtotal++, dock.orient_total++
            )
            {
              if (score.flag)
              {
                dock.score_time -= elapsed_time (NULL);

                if (get_anchor_score
                  (
                    &dock, &label, &score,
                    &mol_ref, &mol_init, &mol_ori, &mol_score,
                    &best_anchors, dock.orient_subtotal,
                    !label.flex.periph_flag && write_flag
                  ) == TRUE)
                {
                  dock.score_subtotal++;
                  dock.score_total++;
                }

                dock.score_time += elapsed_time (NULL);
              }

              else
              {
                if (dock.write_orients)
                {
                  calc_rmsd (&mol_ori, &mol_init);

                  write_ligand
                  (
                    &dock,
                    &score,
                    &mol_ori,
                    score.type[NONE].file_name,
                    score.type[NONE].file
                  );
                }

                update_list (&best_anchors, &mol_ori);
              }
            } /* End of orientation loop */

            dock.orient_time += elapsed_time (NULL);
          }

/*
*         If no screening, scoring, or orienting requested,
*         then write out anchor conformation
*         1/97 te
*/
          else
          {
            write_ligand
            (
              &dock,
              &score,
              &mol_conf,
              score.type[NONE].file_name,
              score.type[NONE].file
            );

            update_list (&best_anchors, &mol_conf);
          }
        }  /* End of anchor conformation loop */

/*
*       From the best scoring anchor configurations, either grow out the
*       remainder of the molecule, or just store the anchors
*       2/97 te
*/
        if (score.flag && label.flex.periph_flag)
        {
          dock.periph_time -= elapsed_time (NULL);

          get_peripheral_score
          (
            &dock, &label, &score,
            &mol_ref, &mol_init, &mol_conf, &mol_ori, &mol_score,
            &best_anchors, &best_orients,
            label.flex.periph_flag && write_flag,
            dock.anchor_total
          );

          dock.periph_time += elapsed_time (NULL);
        }

        else
          merge_lists (&score, &best_orients, &best_anchors);

      }  /* End of anchor loop */

/*
*     Check that an anchor was found
*     6/97 te
*/
      if (dock.anchor_total > 0)
      {
        if (dock.performance_flag && score.minimize.flag)
        {
          score.minimize.call_total += score.minimize.call_sub_total;
          score.minimize.call_max =
            MAX (score.minimize.call_max, score.minimize.call_sub_total);
          score.minimize.call_min =
            MIN (score.minimize.call_min, score.minimize.call_sub_total);
          score.minimize.call_sub_total = 0;
        }

        dock.conform_time += elapsed_time (NULL);
        lig_proc++;

/*
*       Print out the results for ligand
*       11/96 te
*/
        if (score.flag || orient.flag)
          output_score_info
            (&dock, &label, &score, &mol_score, &best_orients, lig_proc - 1,
            elapsed_time (NULL));

        else
        {
          fprintf (global.outfile, ".");

          if (lig_proc % 50 == 0)
            fprintf (global.outfile, "| %d mols, %.4g secs\n",
              lig_proc, elapsed_time (NULL));

          fflush (global.outfile);
        }

/*
*       Write/store this molecule if its orientations weren't written out yet.
*       11/96 te
*/
        if (!write_flag)
        {
          if (score.flag)
          {
/*
*           Either store the best orientation(s)
*           3/96 te
*/
            if (dock.rank_ligands)
            {
              inter_lists (&score, &best_orients, &mol_init);
              merge_lists (&score, &best_ligands, &best_orients);
            }

/*
*           Or write the best orientation(s) out to a file
*           3/96 te
*/
            else
              write_topscorers
              (
                &dock,
                &score,
                &best_orients,
                &mol_ref,
                &mol_out
              );
          }
        }

        else if (label.chemical.screen.construct_flag)
          write_ligand
          (
            &dock,
            &score,