dcdplugin.c

分子动力学计算程序

/***************************************************************************
 *cr                                                                       
 *cr            (C) Copyright 1995-2006 The Board of Trustees of the           
 *cr                        University of Illinois                       
 *cr                         All Rights Reserved                        
 *cr                                                                   
 ***************************************************************************/

/***************************************************************************
 * RCS INFORMATION:
 *
 *      $RCSfile: dcdplugin.c,v $
 *      $Author: jim $       $Locker:  $             $State: Exp $
 *      $Revision: 1.2 $       $Date: 2006/01/20 20:23:55 $
 *
 ***************************************************************************
 * DESCRIPTION:
 *   Code for reading and writing Charmm, NAMD, and X-PLOR format 
 *   molecular dynamic trajectory files.
 *
 * TODO:
 *   Integrate improvements from the NAMD source tree
 *    - NAMD's writer code has better type-correctness for the sizes
 *      of "int".  NAMD uses "int32" explicitly, which is required on
 *      machines like the T3E.  VMD's version of the code doesn't do that
 *      presently.
 *
 *  Try various alternative I/O API options:
 *   - use mmap(), with read-once flags
 *   - use O_DIRECT open mode on new revs of Linux kernel 
 *   - use directio() call on a file descriptor to enable on Solaris
 *   - use aio_open()/read()/write()
 *   - use readv/writev() etc.
 *
 *  Standalone test binary compilation flags:
 *  cc -fast -xarch=v9a -I../../include -DTEST_DCDPLUGIN dcdplugin.c \
 *    -o ~/bin/readdcd -lm
 *
 *  Profiling flags:
 *  cc -xpg -fast -xarch=v9a -g -I../../include -DTEST_DCDPLUGIN dcdplugin.c \
 *    -o ~/bin/readdcd -lm
 *
 ***************************************************************************/

#include "largefiles.h"   /* platform dependent 64-bit file I/O defines */

#include <stdio.h>
#include <sys/stat.h>
#include <sys/types.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include "endianswap.h"
#include "fastio.h"
#include "molfile_plugin.h"

#ifndef M_PI_2
#define M_PI_2 1.57079632679489661922
#endif

typedef struct {
  fio_fd fd;
  int natoms;
  int nsets;
  int setsread;
  int istart;
  int nsavc;
  double delta;
  int nfixed;
  float *x, *y, *z;
  int *freeind;
  float *fixedcoords;
  int reverse;
  int charmm;  
  int first;
  int with_unitcell;
} dcdhandle;

/* Define error codes that may be returned by the DCD routines */
#define DCD_SUCCESS      0  /* No problems                     */
#define DCD_EOF         -1  /* Normal EOF                      */
#define DCD_DNE         -2  /* DCD file does not exist         */
#define DCD_OPENFAILED  -3  /* Open of DCD file failed         */
#define DCD_BADREAD     -4  /* read call on DCD file failed    */
#define DCD_BADEOF      -5  /* premature EOF found in DCD file */
#define DCD_BADFORMAT   -6  /* format of DCD file is wrong     */
#define DCD_FILEEXISTS  -7  /* output file already exists      */
#define DCD_BADMALLOC   -8  /* malloc failed                   */
#define DCD_BADWRITE    -9  /* write call on DCD file failed   */

/* Define feature flags for this DCD file */
#define DCD_IS_CHARMM       0x01
#define DCD_HAS_4DIMS       0x02
#define DCD_HAS_EXTRA_BLOCK 0x04

/* defines used by write_dcdstep */
#define NFILE_POS 8L
#define NSTEP_POS 20L

/* READ Macro to make porting easier */
#define READ(fd, buf, size)  fio_fread(((void *) buf), (size), 1, (fd))

/* WRITE Macro to make porting easier */
#define WRITE(fd, buf, size) fio_fwrite(((void *) buf), (size), 1, (fd))

/* XXX This is broken - fread never returns -1 */
#define CHECK_FREAD(X, msg) if (X==-1) { return(DCD_BADREAD); }
#define CHECK_FEOF(X, msg)  if (X==0)  { return(DCD_BADEOF); }

/*
 * Read the header information from a dcd file.
 * Input: fd - a file struct opened for binary reading.
 * Output: 0 on success, negative error code on failure.
 * Side effects: *natoms set to number of atoms per frame
 *               *nsets set to number of frames in dcd file
 *               *istart set to starting timestep of dcd file
 *               *nsavc set to timesteps between dcd saves
 *               *delta set to value of trajectory timestep
 *               *nfixed set to number of fixed atoms 
 *               *freeind may be set to heap-allocated space
 *               *reverse set to one if reverse-endian, zero if not.
 *               *charmm set to internal code for handling charmm data.
 */
static int read_dcdheader(fio_fd fd, int *N, int *NSET, int *ISTART, 
                   int *NSAVC, double *DELTA, int *NAMNF, 
                   int **FREEINDEXES, float **fixedcoords, int *reverseEndian, 
                   int *charmm)
{
  int input_integer;  /* buffer space */
  int i, ret_val;
  char hdrbuf[84];    /* char buffer used to store header */
  int NTITLE;

  /*  First thing in the file should be an 84 */
  ret_val = READ(fd, &input_integer, sizeof(int));
  CHECK_FREAD(ret_val, "reading first int from dcd file");
  CHECK_FEOF(ret_val, "reading first int from dcd file");

  /* Check magic number in file header and determine byte order*/
  if (input_integer != 84) {
    /* check to see if its merely reversed endianism     */
    /* rather than a totally incorrect file magic number */
    swap4_aligned(&input_integer, 1);

    if (input_integer == 84) {
      *reverseEndian=1;
    } else {
      /* not simply reversed endianism, but something rather more evil */
      return DCD_BADFORMAT;
    }
  } else {
    *reverseEndian=0;    
  }

  /* Buffer the entire header for random access */
  ret_val = READ(fd, hdrbuf, 84);
  CHECK_FREAD(ret_val, "buffering header");
  CHECK_FEOF(ret_val, "buffering header");

  /* Check for the ID string "COORD" */
  if (hdrbuf[0] != 'C' || hdrbuf[1] != 'O' ||
      hdrbuf[2] != 'R' || hdrbuf[3] != 'D') {
    return DCD_BADFORMAT;
  }

  /* CHARMm-genereate DCD files set the last integer in the     */
  /* header, which is unused by X-PLOR, to its version number.  */
  /* Checking if this is nonzero tells us this is a CHARMm file */
  /* and to look for other CHARMm flags.                        */
  if (*((int *) (hdrbuf + 80)) != 0) {
    (*charmm) = DCD_IS_CHARMM;
    if (*((int *) (hdrbuf + 44)) != 0)
      (*charmm) |= DCD_HAS_EXTRA_BLOCK;

    if (*((int *) (hdrbuf + 48)) == 1)
      (*charmm) |= DCD_HAS_4DIMS;
  } else {
    (*charmm) = 0;
  }

  /* Store the number of sets of coordinates (NSET) */
  (*NSET) = *((int *) (hdrbuf + 4));
  if (*reverseEndian) swap4_unaligned(NSET, 1);

  /* Store ISTART, the starting timestep */
  (*ISTART) = *((int *) (hdrbuf + 8));
  if (*reverseEndian) swap4_unaligned(ISTART, 1);

  /* Store NSAVC, the number of timesteps between dcd saves */
  (*NSAVC) = *((int *) (hdrbuf + 12));
  if (*reverseEndian) swap4_unaligned(NSAVC, 1);

  /* Store NAMNF, the number of fixed atoms */
  (*NAMNF) = *((int *) (hdrbuf + 36));
  if (*reverseEndian) swap4_unaligned(NAMNF, 1);

  /* Read in the timestep, DELTA */
  /* Note: DELTA is stored as a double with X-PLOR but as a float with CHARMm */
  if ((*charmm) & DCD_IS_CHARMM) {
    float ftmp;
    ftmp = *((float *)(hdrbuf+40)); /* is this safe on Alpha? */
    if (*reverseEndian)
      swap4_aligned(&ftmp, 1);

    *DELTA = (double)ftmp;
  } else {
    (*DELTA) = *((double *)(hdrbuf + 40));
    if (*reverseEndian) swap8_unaligned(DELTA, 1);
  }

  /* Get the end size of the first block */
  ret_val = READ(fd, &input_integer, sizeof(int));
  CHECK_FREAD(ret_val, "reading second 84 from dcd file");
  CHECK_FEOF(ret_val, "reading second 84 from dcd file");
  if (*reverseEndian) swap4_aligned(&input_integer, 1);

  if (input_integer != 84) {
    return DCD_BADFORMAT;
  }

  /* Read in the size of the next block */
  ret_val = READ(fd, &input_integer, sizeof(int));
  CHECK_FREAD(ret_val, "reading size of title block");
  CHECK_FEOF(ret_val, "reading size of title block");
  if (*reverseEndian) swap4_aligned(&input_integer, 1);

  if (((input_integer-4) % 80) == 0) {
    /* Read NTITLE, the number of 80 character title strings there are */
    ret_val = READ(fd, &NTITLE, sizeof(int));
    CHECK_FREAD(ret_val, "reading NTITLE");
    CHECK_FEOF(ret_val, "reading NTITLE");
    if (*reverseEndian) swap4_aligned(&NTITLE, 1);

    for (i=0; i<NTITLE; i++) {
      fio_fseek(fd, 80, FIO_SEEK_CUR);
      CHECK_FEOF(ret_val, "reading TITLE");
    }

    /* Get the ending size for this block */
    ret_val = READ(fd, &input_integer, sizeof(int));

    CHECK_FREAD(ret_val, "reading size of title block");
    CHECK_FEOF(ret_val, "reading size of title block");
  } else {
    return DCD_BADFORMAT;
  }

  /* Read in an integer '4' */
  ret_val = READ(fd, &input_integer, sizeof(int));
  CHECK_FREAD(ret_val, "reading a '4'");
  CHECK_FEOF(ret_val, "reading a '4'");
  if (*reverseEndian) swap4_aligned(&input_integer, 1);

  if (input_integer != 4) {
    return DCD_BADFORMAT;
  }

  /* Read in the number of atoms */
  ret_val = READ(fd, N, sizeof(int));
  CHECK_FREAD(ret_val, "reading number of atoms");
  CHECK_FEOF(ret_val, "reading number of atoms");
  if (*reverseEndian) swap4_aligned(N, 1);

  /* Read in an integer '4' */
  ret_val = READ(fd, &input_integer, sizeof(int));
  CHECK_FREAD(ret_val, "reading a '4'");
  CHECK_FEOF(ret_val, "reading a '4'");
  if (*reverseEndian) swap4_aligned(&input_integer, 1);

  if (input_integer != 4) {
    return DCD_BADFORMAT;
  }

  *FREEINDEXES = NULL;
  *fixedcoords = NULL;
  if (*NAMNF != 0) {
    (*FREEINDEXES) = (int *) calloc(((*N)-(*NAMNF)), sizeof(int));
    if (*FREEINDEXES == NULL)
      return DCD_BADMALLOC;

    *fixedcoords = (float *) calloc((*N)*4 - (*NAMNF), sizeof(float));
    if (*fixedcoords == NULL)
      return DCD_BADMALLOC;

    /* Read in index array size */
    ret_val = READ(fd, &input_integer, sizeof(int));
    CHECK_FREAD(ret_val, "reading size of index array");
    CHECK_FEOF(ret_val, "reading size of index array");
    if (*reverseEndian) swap4_aligned(&input_integer, 1);

    if (input_integer != ((*N)-(*NAMNF))*4) {
      return DCD_BADFORMAT;
    }

    ret_val = READ(fd, (*FREEINDEXES), ((*N)-(*NAMNF))*sizeof(int));
    CHECK_FREAD(ret_val, "reading size of index array");
    CHECK_FEOF(ret_val, "reading size of index array");

    if (*reverseEndian)
      swap4_aligned((*FREEINDEXES), ((*N)-(*NAMNF)));

    ret_val = READ(fd, &input_integer, sizeof(int));
    CHECK_FREAD(ret_val, "reading size of index array");
    CHECK_FEOF(ret_val, "reading size of index array");
    if (*reverseEndian) swap4_aligned(&input_integer, 1);

    if (input_integer != ((*N)-(*NAMNF))*4) {
      return DCD_BADFORMAT;
    }
  }

  return DCD_SUCCESS;
}

static int read_charmm_extrablock(fio_fd fd, int charmm, int reverseEndian,
                                  float *unitcell) {
  int i, input_integer;

  if ((charmm & DCD_IS_CHARMM) && (charmm & DCD_HAS_EXTRA_BLOCK)) {
    /* Leading integer must be 48 */
    if (fio_fread(&input_integer, sizeof(int), 1, fd) != 1)
      return DCD_BADREAD; 
    if (reverseEndian) swap4_aligned(&input_integer, 1);
    if (input_integer == 48) {
      double tmp[6];
      if (fio_fread(tmp, 48, 1, fd) != 1) return DCD_BADREAD;
      if (reverseEndian) 
        swap8_aligned(tmp, 6);
      for (i=0; i<6; i++) unitcell[i] = (float)tmp[i];
    } else {
      /* unrecognized block, just skip it */
      if (fio_fseek(fd, input_integer, FIO_SEEK_CUR)) return DCD_BADREAD;
    }
    if (fio_fread(&input_integer, sizeof(int), 1, fd) != 1) return DCD_BADREAD; 
  } 

  return DCD_SUCCESS;
}

static int read_fixed_atoms(fio_fd fd, int N, int num_free, const int *indexes,
                            int reverseEndian, const float *fixedcoords, 
                            float *freeatoms, float *pos) {
  int i, input_integer;
  
  /* Read leading integer */
  if (fio_fread(&input_integer, sizeof(int), 1, fd) != 1) return DCD_BADREAD;
  if (reverseEndian) swap4_aligned(&input_integer, 1);
  if (input_integer != 4*num_free) return DCD_BADFORMAT;