dcdplugin.c

分子动力学计算程序

  
  /* Read free atom coordinates */
  if (fio_fread(freeatoms, 4*num_free, 1, fd) != 1) return DCD_BADREAD;
  if (reverseEndian)
    swap4_aligned(freeatoms, num_free);

  /* Copy fixed and free atom coordinates into position buffer */
  memcpy(pos, fixedcoords, 4*N);
  for (i=0; i<num_free; i++)
    pos[indexes[i]-1] = freeatoms[i];

  /* Read trailing 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;

  return DCD_SUCCESS;
}
  
static int read_charmm_4dim(fio_fd fd, int charmm, int reverseEndian) {
  int input_integer;

  /* If this is a CHARMm file and contains a 4th dimension block, */
  /* we must skip past it to avoid problems                       */
  if ((charmm & DCD_IS_CHARMM) && (charmm & DCD_HAS_4DIMS)) {
    if (fio_fread(&input_integer, sizeof(int), 1, fd) != 1) return DCD_BADREAD;  
    if (reverseEndian) swap4_aligned(&input_integer, 1);
    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;
}

/* 
 * Read a dcd timestep from a dcd file
 * Input: fd - a file struct opened for binary reading, from which the 
 *             header information has already been read.
 *        natoms, nfixed, first, *freeind, reverse, charmm - the corresponding 
 *             items as set by read_dcdheader
 *        first - true if this is the first frame we are reading.
 *        x, y, z: space for natoms each of floats.
 *        unitcell - space for six floats to hold the unit cell data.  
 *                   Not set if no unit cell data is present.
 * Output: 0 on success, negative error code on failure.
 * Side effects: x, y, z contain the coordinates for the timestep read.
 *               unitcell holds unit cell data if present.
 */
static int read_dcdstep(fio_fd fd, int N, float *X, float *Y, float *Z, 
                        float *unitcell, int num_fixed,
                        int first, int *indexes, float *fixedcoords, 
                        int reverseEndian, int charmm) {
  int ret_val;   /* Return value from read */

  if ((num_fixed==0) || first) {
    int tmpbuf[6];      /* temp storage for reading formatting info */
    fio_iovec iov[7];   /* I/O vector for fio_readv() call          */
    fio_size_t readlen; /* number of bytes actually read            */
    int i;

    /* if there are no fixed atoms or this is the first timestep read */
    /* then we read all coordinates normally.                         */

    /* read the charmm periodic cell information */
    /* XXX this too should be read together with the other items in a */
    /*     single fio_readv() call in order to prevent lots of extra  */
    /*     kernel/user context switches.                              */
    ret_val = read_charmm_extrablock(fd, charmm, reverseEndian, unitcell);
    if (ret_val) return ret_val;

    /* setup the I/O vector for the call to fio_readv() */
    iov[0].iov_base = (fio_caddr_t) &tmpbuf[0]; /* read format integer    */
    iov[0].iov_len  = sizeof(int);

    iov[1].iov_base = (fio_caddr_t) X;          /* read X coordinates     */
    iov[1].iov_len  = sizeof(float)*N;

    iov[2].iov_base = (fio_caddr_t) &tmpbuf[1]; /* read 2 format integers */
    iov[2].iov_len  = sizeof(int) * 2;

    iov[3].iov_base = (fio_caddr_t) Y;          /* read Y coordinates     */
    iov[3].iov_len  = sizeof(float)*N;

    iov[4].iov_base = (fio_caddr_t) &tmpbuf[3]; /* read 2 format integers */
    iov[4].iov_len  = sizeof(int) * 2;

    iov[5].iov_base = (fio_caddr_t) Z;          /* read Y coordinates     */
    iov[5].iov_len  = sizeof(float)*N;

    iov[6].iov_base = (fio_caddr_t) &tmpbuf[5]; /* read format integer    */
    iov[6].iov_len  = sizeof(int);

    readlen = fio_readv(fd, &iov[0], 7);

    if (readlen != (6*sizeof(int) + 3*N*sizeof(float)))
      return DCD_BADREAD;

    /* convert endianism if necessary */
    if (reverseEndian) {
      swap4_aligned(&tmpbuf[0], 6);
      swap4_aligned(X, N);
      swap4_aligned(Y, N);
      swap4_aligned(Z, N);
    }

    /* double-check the fortran format size values for safety */
    for (i=0; i<6; i++) {
      if (tmpbuf[i] != sizeof(float)*N) return DCD_BADFORMAT;
    }

    /* copy fixed atom coordinates into fixedcoords array if this was the */
    /* first timestep, to be used from now on.  We just copy all atoms.   */
    if (num_fixed && first) {
      memcpy(fixedcoords, X, N*sizeof(float));
      memcpy(fixedcoords+N, Y, N*sizeof(float));
      memcpy(fixedcoords+2*N, Z, N*sizeof(float));
    }

    /* read in the optional charmm 4th array */
    /* XXX this too should be read together with the other items in a */
    /*     single fio_readv() call in order to prevent lots of extra  */
    /*     kernel/user context switches.                              */
    ret_val = read_charmm_4dim(fd, charmm, reverseEndian);
    if (ret_val) return ret_val;
  } else {
    /* if there are fixed atoms, and this isn't the first frame, then we */
    /* only read in the non-fixed atoms for all subsequent timesteps.    */
    ret_val = read_charmm_extrablock(fd, charmm, reverseEndian, unitcell);
    if (ret_val) return ret_val;
    ret_val = read_fixed_atoms(fd, N, N-num_fixed, indexes, reverseEndian,
                               fixedcoords, fixedcoords+3*N, X);
    if (ret_val) return ret_val;
    ret_val = read_fixed_atoms(fd, N, N-num_fixed, indexes, reverseEndian,
                               fixedcoords+N, fixedcoords+3*N, Y);
    if (ret_val) return ret_val;
    ret_val = read_fixed_atoms(fd, N, N-num_fixed, indexes, reverseEndian,
                               fixedcoords+2*N, fixedcoords+3*N, Z);
    if (ret_val) return ret_val;
    ret_val = read_charmm_4dim(fd, charmm, reverseEndian);
    if (ret_val) return ret_val;
  }

  return DCD_SUCCESS;
}


/* 
 * Skip past a timestep.  If there are fixed atoms, this cannot be used with
 * the first timestep.  
 * Input: fd - a file struct from which the header has already been read
 *        natoms - number of atoms per timestep
 *        nfixed - number of fixed atoms
 *        charmm - charmm flags as returned by read_dcdheader
 * Output: 0 on success, negative error code on failure.
 * Side effects: One timestep will be skipped; fd will be positioned at the
 *               next timestep.
 */
static int skip_dcdstep(fio_fd fd, int natoms, int nfixed, int charmm) {
  int seekoffset = 0;

  /* Skip charmm extra block */
  if ((charmm & DCD_IS_CHARMM) && (charmm & DCD_HAS_EXTRA_BLOCK)) {
    seekoffset += 4 + 48 + 4;
  }

  /* For each atom set, seek past an int, the free atoms, and another int. */
  seekoffset += 3 * (2 + natoms - nfixed) * 4;

  /* Assume that charmm 4th dim is the same size as the other three. */
  if ((charmm & DCD_IS_CHARMM) && (charmm & DCD_HAS_4DIMS)) {
    seekoffset += (2 + natoms - nfixed) * 4;
  }
 
  if (fio_fseek(fd, seekoffset, FIO_SEEK_CUR)) return DCD_BADEOF;

  return DCD_SUCCESS;
}


/* 
 * Write a timestep to a dcd file
 * Input: fd - a file struct for which a dcd header has already been written
 *       curframe: Count of frames written to this file, starting with 1.
 *       curstep: Count of timesteps elapsed = istart + curframe * nsavc.
 *        natoms - number of elements in x, y, z arrays
 *        x, y, z: pointers to atom coordinates
 * Output: 0 on success, negative error code on failure.
 * Side effects: coordinates are written to the dcd file.
 */
static int write_dcdstep(fio_fd fd, int curframe, int curstep, int N, 
                  const float *X, const float *Y, const float *Z, 
                  const double *unitcell, int charmm) {
  int out_integer;
  int rc;

  if (charmm) {
    /* write out optional unit cell */
    if (unitcell != NULL) {
      out_integer = 48; /* 48 bytes (6 doubles) */
      fio_write_int32(fd, out_integer);
      WRITE(fd, unitcell, out_integer);
      fio_write_int32(fd, out_integer);
    }
  }

  /* write out coordinates */
  out_integer = N*4; /* N*4 bytes per X/Y/Z array (N floats per array) */
  fio_write_int32(fd, out_integer);
  if (fio_fwrite((void *) X, out_integer, 1, fd) != 1) return DCD_BADWRITE;
  fio_write_int32(fd, out_integer);
  fio_write_int32(fd, out_integer);
  if (fio_fwrite((void *) Y, out_integer, 1, fd) != 1) return DCD_BADWRITE;
  fio_write_int32(fd, out_integer);
  fio_write_int32(fd, out_integer);
  if (fio_fwrite((void *) Z, out_integer, 1, fd) != 1) return DCD_BADWRITE;
  fio_write_int32(fd, out_integer);

  /* update the DCD header information */
  fio_fseek(fd, NFILE_POS, FIO_SEEK_SET);
  fio_write_int32(fd, curframe);
  fio_fseek(fd, NSTEP_POS, FIO_SEEK_SET);
  fio_write_int32(fd, curstep);
  fio_fseek(fd, 0, FIO_SEEK_END);

  return DCD_SUCCESS;
}

/*
 * Write a header for a new dcd file
 * Input: fd - file struct opened for binary writing
 *        remarks - string to be put in the remarks section of the header.  
 *                  The string will be truncated to 70 characters.
 *        natoms, istart, nsavc, delta - see comments in read_dcdheader
 * Output: 0 on success, negative error code on failure.
 * Side effects: Header information is written to the dcd file.
 */
static int write_dcdheader(fio_fd fd, const char *remarks, int N, 
                    int ISTART, int NSAVC, double DELTA, int with_unitcell,
                    int charmm) {
  int out_integer;
  float out_float;
  char title_string[200];
  time_t cur_time;
  struct tm *tmbuf;
  char time_str[81];

  out_integer = 84;
  WRITE(fd, (char *) & out_integer, sizeof(int));
  strcpy(title_string, "CORD");
  WRITE(fd, title_string, 4);
  fio_write_int32(fd, 0);      /* Number of frames in file, none written yet   */
  fio_write_int32(fd, ISTART); /* Starting timestep                            */
  fio_write_int32(fd, NSAVC);  /* Timesteps between frames written to the file */
  fio_write_int32(fd, 0);      /* Number of timesteps in simulation            */
  fio_write_int32(fd, 0);      /* NAMD writes NSTEP or ISTART - NSAVC here?    */
  fio_write_int32(fd, 0);
  fio_write_int32(fd, 0);
  fio_write_int32(fd, 0);
  fio_write_int32(fd, 0);
  if (charmm) {
    out_float = DELTA;
    WRITE(fd, (char *) &out_float, sizeof(float));
    if (with_unitcell) {
      fio_write_int32(fd, 1);
    } else {
      fio_write_int32(fd, 0);
    }
  } else {
    WRITE(fd, (char *) &DELTA, sizeof(double));
  }
  fio_write_int32(fd, 0);
  fio_write_int32(fd, 0);
  fio_write_int32(fd, 0);
  fio_write_int32(fd, 0);
  fio_write_int32(fd, 0);
  fio_write_int32(fd, 0);
  fio_write_int32(fd, 0);
  fio_write_int32(fd, 0);
  if (charmm) {
    fio_write_int32(fd, 24); /* Pretend to be Charmm version 24 */
  } else {
    fio_write_int32(fd, 0);
  }
  fio_write_int32(fd, 84);
  fio_write_int32(fd, 164);
  fio_write_int32(fd, 2);

  strncpy(title_string, remarks, 80);
  title_string[79] = '\0';
  WRITE(fd, title_string, 80);

  cur_time=time(NULL);
  tmbuf=localtime(&cur_time);
  strftime(time_str, 80, "REMARKS Created %d %B, %Y at %R", tmbuf);
  WRITE(fd, time_str, 80);

  fio_write_int32(fd, 164);
  fio_write_int32(fd, 4);
  fio_write_int32(fd, N);
  fio_write_int32(fd, 4);

  return DCD_SUCCESS;
}


/*
 * clean up dcd data
 * Input: nfixed, freeind - elements as returned by read_dcdheader
 * Output: None
 * Side effects: Space pointed to by freeind is freed if necessary.
 */
static void close_dcd_read(int *indexes, float *fixedcoords) {
  free(indexes);
  free(fixedcoords);
}





static void *open_dcd_read(const char *path, const char *filetype, 
    int *natoms) {
  dcdhandle *dcd;
  fio_fd fd;
  int rc;
  struct stat stbuf;

  if (!path) return NULL;

  /* See if the file exists, and get its size */
  memset(&stbuf, 0, sizeof(struct stat));
  if (stat(path, &stbuf)) {
    fprintf(stderr, "Could not access file '%s'.\n", path);
    return NULL;
  }

  if (fio_open(path, FIO_READ, &fd) < 0) {
    fprintf(stderr, "Could not open file '%s' for reading.\n", path);
    return NULL;
  }

  dcd = (dcdhandle *)malloc(sizeof(dcdhandle));
  memset(dcd, 0, sizeof(dcdhandle));
  dcd->fd = fd;

  if ((rc = read_dcdheader(dcd->fd, &dcd->natoms, &dcd->nsets, &dcd->istart, 
         &dcd->nsavc, &dcd->delta, &dcd->nfixed, &dcd->freeind, 
         &dcd->fixedcoords, &dcd->reverse, &dcd->charmm))) {
    fprintf(stderr, "read_dcdheader returned %d\n", rc);