Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Related Pages  

dcdplugin.c

Go to the documentation of this file.
00001 /***************************************************************************
00002  *cr                                                                       
00003  *cr            (C) Copyright 1995-2016 The Board of Trustees of the           
00004  *cr                        University of Illinois                       
00005  *cr                         All Rights Reserved                        
00006  *cr                                                                   
00007  ***************************************************************************/
00008 
00009 /***************************************************************************
00010  * RCS INFORMATION:
00011  *
00012  *      $RCSfile: dcdplugin.c,v $
00013  *      $Author: johns $       $Locker:  $             $State: Exp $
00014  *      $Revision: 1.88 $       $Date: 2020/12/17 17:14:07 $
00015  *
00016  ***************************************************************************
00017  * DESCRIPTION:
00018  *   Code for reading and writing CHARMM, NAMD, and X-PLOR format 
00019  *   molecular dynamic trajectory files.
00020  *
00021  * TODO:
00022  *   Integrate improvements from the NAMD source tree
00023  *    - NAMD's writer code has better type-correctness for the sizes
00024  *      of "int".  NAMD uses "int32" explicitly, which is required on
00025  *      machines like the T3E.  VMD's version of the code doesn't do that
00026  *      presently.
00027  *
00028  *  Try various alternative I/O API options:
00029  *   - use mmap(), with read-once flags
00030  *   - use O_DIRECT open mode on new revs of Linux kernel 
00031  *   - use directio() call on a file descriptor to enable on Solaris
00032  *   - use aio_open()/read()/write()
00033  *   - use readv/writev() etc.
00034  *
00035  *  Standalone test binary compilation flags:
00036  *  cc -fast -xarch=v9a -I../../include -DTEST_DCDPLUGIN dcdplugin.c \
00037  *    -o ~/bin/readdcd -lm
00038  *
00039  *  Profiling flags:
00040  *  cc -xpg -fast -xarch=v9a -g -I../../include -DTEST_DCDPLUGIN dcdplugin.c \
00041  *    -o ~/bin/readdcd -lm
00042  *
00043  ***************************************************************************/
00044 
00045 #include "largefiles.h"   /* platform dependent 64-bit file I/O defines */
00046 #include "fastio.h"       /* must come before others, for O_DIRECT...   */
00047 
00048 #include <stdio.h>
00049 #include <sys/stat.h>
00050 #include <sys/types.h>
00051 #include <stdlib.h>
00052 #include <stddef.h>
00053 #include <string.h>
00054 #include <math.h>
00055 #include <time.h>
00056 #include "endianswap.h"
00057 #include "molfile_plugin.h"
00058 
00059 #ifndef M_PI_2
00060 #define M_PI_2 1.57079632679489661922
00061 #endif
00062 
00063 #define RECSCALE32BIT 1
00064 #define RECSCALE64BIT 2
00065 #define RECSCALEMAX   2
00066 
00067 typedef struct {
00068   fio_fd fd;
00069   int natoms;
00070   int nsets;
00071   int setsread;
00072   int istart;
00073   int nsavc;
00074   double delta;
00075   int nfixed;
00076   float *x, *y, *z;
00077   int *freeind;
00078   float *fixedcoords;
00079   int reverse;
00080   int charmm;  
00081   int first;
00082   int with_unitcell;
00083 } dcdhandle;
00084 
00085 /* Define error codes that may be returned by the DCD routines */
00086 #define DCD_SUCCESS      0  /* No problems                     */
00087 #define DCD_EOF         -1  /* Normal EOF                      */
00088 #define DCD_DNE         -2  /* DCD file does not exist         */
00089 #define DCD_OPENFAILED  -3  /* Open of DCD file failed         */
00090 #define DCD_BADREAD     -4  /* read call on DCD file failed    */
00091 #define DCD_BADEOF      -5  /* premature EOF found in DCD file */
00092 #define DCD_BADFORMAT   -6  /* format of DCD file is wrong     */
00093 #define DCD_FILEEXISTS  -7  /* output file already exists      */
00094 #define DCD_BADMALLOC   -8  /* malloc failed                   */
00095 #define DCD_BADWRITE    -9  /* write call on DCD file failed   */
00096 
00097 /* Define feature flags for this DCD file */
00098 #define DCD_IS_XPLOR        0x00
00099 #define DCD_IS_CHARMM       0x01
00100 #define DCD_HAS_4DIMS       0x02
00101 #define DCD_HAS_EXTRA_BLOCK 0x04
00102 #define DCD_HAS_64BIT_REC   0x08
00103 
00104 /* defines used by write_dcdstep */
00105 #define NFILE_POS 8L
00106 #define NSTEP_POS 20L
00107 
00108 /* READ Macro to make porting easier */
00109 #define READ(fd, buf, size)  fio_fread(((void *) buf), (size), 1, (fd))
00110 
00111 /* WRITE Macro to make porting easier */
00112 #define WRITE(fd, buf, size) fio_fwrite(((void *) buf), (size), 1, (fd))
00113 
00114 /* XXX This is broken - fread never returns -1 */
00115 #define CHECK_FREAD(X, msg) if (X==-1) { return(DCD_BADREAD); }
00116 #define CHECK_FEOF(X, msg)  if (X==0)  { return(DCD_BADEOF); }
00117 
00118 
00119 /* print DCD error in a human readable way */
00120 static void print_dcderror(const char *func, int errcode) {
00121   const char *errstr;
00122 
00123   switch (errcode) {
00124     case DCD_EOF:         errstr = "end of file"; break;
00125     case DCD_DNE:         errstr = "file not found"; break;
00126     case DCD_OPENFAILED:  errstr = "file open failed"; break;
00127     case DCD_BADREAD:     errstr = "error during read"; break;
00128     case DCD_BADEOF:      errstr = "premature end of file"; break;
00129     case DCD_BADFORMAT:   errstr = "corruption or unrecognized file structure"; break;
00130     case DCD_FILEEXISTS:  errstr = "output file already exists"; break;
00131     case DCD_BADMALLOC:   errstr = "memory allocation failed"; break;
00132     case DCD_BADWRITE:    errstr = "error during write"; break;
00133     case DCD_SUCCESS:     
00134     default:
00135       errstr = "no error";
00136       break;
00137   } 
00138   printf("dcdplugin) %s: %s\n", func, errstr); 
00139 }
00140 
00141 
00142 /*
00143  * Read the header information from a dcd file.
00144  * Input: fd - a file struct opened for binary reading.
00145  * Output: 0 on success, negative error code on failure.
00146  * Side effects: *natoms set to number of atoms per frame
00147  *               *nsets set to number of frames in dcd file
00148  *               *istart set to starting timestep of dcd file
00149  *               *nsavc set to timesteps between dcd saves
00150  *               *delta set to value of trajectory timestep
00151  *               *nfixed set to number of fixed atoms 
00152  *               *freeind may be set to heap-allocated space
00153  *               *reverse set to one if reverse-endian, zero if not.
00154  *               *charmm set to internal code for handling charmm data.
00155  */
00156 static int read_dcdheader(fio_fd fd, int *N, int *NSET, int *ISTART, 
00157                    int *NSAVC, double *DELTA, int *NAMNF, 
00158                    int **FREEINDEXES, float **fixedcoords, int *reverseEndian, 
00159                    int *charmm)
00160 {
00161   unsigned int input_integer[2];  /* buffer space */
00162   int i, ret_val, rec_scale;
00163   char hdrbuf[84];    /* char buffer used to store header */
00164   int NTITLE;
00165   int dcdcordmagic;
00166   int hugefile = 0;
00167   char *corp = (char *) &dcdcordmagic;
00168 
00169   /* coordinate dcd file magic string 'CORD' */
00170   corp[0] = 'C';
00171   corp[1] = 'O';
00172   corp[2] = 'R';
00173   corp[3] = 'D';
00174 
00175   /* First thing in the file should be an 84.
00176    * some 64-bit compiles have a 64-bit record length indicator,
00177    * so we have to read two ints and check in a more complicated 
00178    * way. :-( */
00179   ret_val = READ(fd, input_integer, 2*sizeof(unsigned int));
00180   CHECK_FREAD(ret_val, "reading first int from dcd file");
00181   CHECK_FEOF(ret_val, "reading first int from dcd file");
00182 
00183   /* Check magic number in file header and determine byte order*/
00184   if ((input_integer[0]+input_integer[1]) == 84) {
00185     *reverseEndian=0;
00186     rec_scale=RECSCALE64BIT;
00187     printf("dcdplugin) detected CHARMM -i8 64-bit DCD file of native endianness\n");
00188   } else if (input_integer[0] == 84 && input_integer[1] == dcdcordmagic) {
00189     *reverseEndian=0;
00190     rec_scale=RECSCALE32BIT;
00191     printf("dcdplugin) detected standard 32-bit DCD file of native endianness\n");
00192   } else {
00193     /* now try reverse endian */
00194     swap4_aligned(input_integer, 2); /* will have to unswap magic if 32-bit */
00195     if ((input_integer[0]+input_integer[1]) == 84) {
00196       *reverseEndian=1;
00197       rec_scale=RECSCALE64BIT;
00198       printf("dcdplugin) detected CHARMM -i8 64-bit DCD file of opposite endianness\n");
00199     } else {
00200       swap4_aligned(&input_integer[1], 1); /* unswap magic (see above) */
00201       if (input_integer[0] == 84 && input_integer[1] == dcdcordmagic) {
00202         *reverseEndian=1;
00203         rec_scale=RECSCALE32BIT;
00204         printf("dcdplugin) detected standard 32-bit DCD file of opposite endianness\n");
00205       } else {
00206         /* not simply reversed endianism or -i8, something rather more evil */
00207         printf("dcdplugin) unrecognized DCD header:\n");
00208         printf("dcdplugin)   [0]: %10d  [1]: %10d\n", input_integer[0], input_integer[1]);
00209         printf("dcdplugin)   [0]: 0x%08x  [1]: 0x%08x\n", input_integer[0], input_integer[1]);
00210         return DCD_BADFORMAT;
00211 
00212       }
00213     }
00214   }
00215 
00216   /* check for magic string, in case of long record markers */
00217   if (rec_scale == RECSCALE64BIT) { 
00218     ret_val = READ(fd, input_integer, sizeof(unsigned int));
00219     if (input_integer[0] != dcdcordmagic) {
00220       printf("dcdplugin) failed to find CORD magic in CHARMM -i8 64-bit DCD file\n");
00221       return DCD_BADFORMAT;
00222     }
00223   }
00224 
00225   /* Buffer the entire header for random access */
00226   ret_val = READ(fd, hdrbuf, 80);
00227   CHECK_FREAD(ret_val, "buffering header");
00228   CHECK_FEOF(ret_val, "buffering header");
00229 
00230   /* CHARMm-genereate DCD files set the last integer in the     */
00231   /* header, which is unused by X-PLOR, to its version number.  */
00232   /* Checking if this is nonzero tells us this is a CHARMm file */
00233   /* and to look for other CHARMm flags.                        */
00234   if (*((int *) (hdrbuf + 76)) != 0) {
00235     (*charmm) = DCD_IS_CHARMM;
00236     if (*((int *) (hdrbuf + 40)) != 0)
00237       (*charmm) |= DCD_HAS_EXTRA_BLOCK;
00238 
00239     if (*((int *) (hdrbuf + 44)) == 1)
00240       (*charmm) |= DCD_HAS_4DIMS;
00241 
00242     if (rec_scale == RECSCALE64BIT)
00243       (*charmm) |= DCD_HAS_64BIT_REC;
00244   
00245   } else {
00246     (*charmm) = DCD_IS_XPLOR; /* must be an X-PLOR format DCD file */
00247   }
00248 
00249   if (*charmm & DCD_IS_CHARMM) {
00250     /* CHARMM and NAMD versions 2.1b1 and later */
00251     printf("dcdplugin) CHARMM format DCD file (also NAMD 2.1 and later)\n");
00252   } else {
00253     /* CHARMM and NAMD versions prior to 2.1b1  */
00254     printf("dcdplugin) X-PLOR format DCD file (also NAMD 2.0 and earlier)\n");
00255   }
00256 
00257   /* Store the number of sets of coordinates (NSET) */
00258   (*NSET) = *((int *) (hdrbuf));
00259   if (*reverseEndian) swap4_unaligned(NSET, 1);
00260 
00261   /* Store ISTART, the starting timestep */
00262   (*ISTART) = *((int *) (hdrbuf + 4));
00263   if (*reverseEndian) swap4_unaligned(ISTART, 1);
00264 
00265   /* Store NSAVC, the number of timesteps between dcd saves */
00266   (*NSAVC) = *((int *) (hdrbuf + 8));
00267   if (*reverseEndian) swap4_unaligned(NSAVC, 1);
00268 
00269   /* Store NAMNF, the number of fixed atoms */
00270   (*NAMNF) = *((int *) (hdrbuf + 32));
00271   if (*reverseEndian) swap4_unaligned(NAMNF, 1);
00272 
00273   /* Read in the timestep, DELTA */
00274   /* Note: DELTA is stored as a double with X-PLOR but as a float with CHARMm */
00275   if ((*charmm) & DCD_IS_CHARMM) {
00276     float ftmp;
00277     ftmp = *((float *)(hdrbuf+36)); /* is this safe on Alpha? */
00278     if (*reverseEndian)
00279       swap4_aligned(&ftmp, 1);
00280 
00281     *DELTA = (double)ftmp;
00282   } else {
00283     (*DELTA) = *((double *)(hdrbuf + 36));
00284     if (*reverseEndian) swap8_unaligned(DELTA, 1);
00285   }
00286 
00287   /* Get the end size of the first block */
00288   ret_val = READ(fd, input_integer, rec_scale*sizeof(int));
00289   CHECK_FREAD(ret_val, "reading second 84 from dcd file");
00290   CHECK_FEOF(ret_val, "reading second 84 from dcd file");
00291   if (*reverseEndian) swap4_aligned(input_integer, rec_scale);
00292 
00293   if (rec_scale == RECSCALE64BIT) {
00294     if ((input_integer[0]+input_integer[1]) != 84) {
00295       return DCD_BADFORMAT;
00296     }
00297   } else {
00298     if (input_integer[0] != 84) {
00299       return DCD_BADFORMAT;
00300     }
00301   }
00302   
00303   /* Read in the size of the next block */
00304   input_integer[1] = 0;
00305   ret_val = READ(fd, input_integer, rec_scale*sizeof(int));
00306   CHECK_FREAD(ret_val, "reading size of title block");
00307   CHECK_FEOF(ret_val, "reading size of title block");
00308   if (*reverseEndian) swap4_aligned(input_integer, rec_scale);
00309 
00310   if ((((input_integer[0]+input_integer[1])-4) % 80) == 0) {
00311     /* Read NTITLE, the number of 80 character title strings there are */
00312     ret_val = READ(fd, &NTITLE, sizeof(int));
00313     CHECK_FREAD(ret_val, "reading NTITLE");
00314     CHECK_FEOF(ret_val, "reading NTITLE");
00315     if (*reverseEndian) swap4_aligned(&NTITLE, 1);
00316 
00317     if (NTITLE < 0) {
00318       printf("dcdplugin) WARNING: Bogus NTITLE value: %d (hex: %08x)\n", 
00319              NTITLE, NTITLE);
00320       return DCD_BADFORMAT;
00321     }
00322 
00323     if (NTITLE > 1000) {
00324       printf("dcdplugin) WARNING: Bogus NTITLE value: %d (hex: %08x)\n", 
00325              NTITLE, NTITLE);
00326       if (NTITLE == 1095062083) {
00327         printf("dcdplugin) WARNING: Broken Vega ZZ 2.4.0 DCD file detected\n");
00328         printf("dcdplugin) Assuming 2 title lines, good luck...\n");
00329         NTITLE = 2;
00330       } else {
00331         printf("dcdplugin) Assuming zero title lines, good luck...\n");
00332         NTITLE = 0;
00333       }
00334     }
00335 
00336     for (i=0; i<NTITLE; i++) {
00337       fio_fseek(fd, 80, FIO_SEEK_CUR);
00338       CHECK_FEOF(ret_val, "reading TITLE");
00339     }
00340 
00341     /* Get the ending size for this block */
00342     ret_val = READ(fd, input_integer, rec_scale*sizeof(int));
00343     CHECK_FREAD(ret_val, "reading size of title block");
00344     CHECK_FEOF(ret_val, "reading size of title block");
00345   } else {
00346     return DCD_BADFORMAT;
00347   }
00348 
00349   /* Read in an integer '4' */
00350   input_integer[1] = 0;
00351   ret_val = READ(fd, input_integer, rec_scale*sizeof(int));
00352   
00353   CHECK_FREAD(ret_val, "reading a '4'");
00354   CHECK_FEOF(ret_val, "reading a '4'");
00355   if (*reverseEndian) swap4_aligned(input_integer, rec_scale);
00356 
00357   if ((input_integer[0]+input_integer[1]) != 4) {
00358     return DCD_BADFORMAT;
00359   }
00360 
00361   /* Read in the number of atoms */
00362   ret_val = READ(fd, N, sizeof(int));
00363   CHECK_FREAD(ret_val, "reading number of atoms");
00364   CHECK_FEOF(ret_val, "reading number of atoms");
00365   if (*reverseEndian) swap4_aligned(N, 1);
00366 
00367   if (*N > (1L<<30)) {
00368     hugefile=1;
00369     printf("dcdplugin) ***\n");
00370     printf("dcdplugin) *** Trajectory contains over 2^30 atoms.\n");
00371     printf("dcdplugin) *** Huge file integer wraparound handling enabled.\n");
00372     printf("dcdplugin) ***\n");
00373   }
00374 
00375   /* Read in an integer '4' */
00376   input_integer[1] = 0;
00377   ret_val = READ(fd, input_integer, rec_scale*sizeof(int));
00378   CHECK_FREAD(ret_val, "reading a '4'");
00379   CHECK_FEOF(ret_val, "reading a '4'");
00380   if (*reverseEndian) swap4_aligned(input_integer, rec_scale);
00381 
00382   if ((input_integer[0]+input_integer[1]) != 4) {
00383     return DCD_BADFORMAT;
00384   }
00385 
00386   *FREEINDEXES = NULL;
00387   *fixedcoords = NULL;
00388   if (*NAMNF != 0) {
00389     (*FREEINDEXES) = (int *) calloc((((ptrdiff_t)(*N))-(*NAMNF)), sizeof(int));
00390     if (*FREEINDEXES == NULL)
00391       return DCD_BADMALLOC;
00392 
00393     *fixedcoords = (float *) calloc(((ptrdiff_t)(*N))*4L - (*NAMNF), sizeof(float));
00394     if (*fixedcoords == NULL)
00395       return DCD_BADMALLOC;
00396 
00397     /* Read in index array size */
00398     input_integer[1]=0;
00399     ret_val = READ(fd, input_integer, rec_scale*sizeof(int));
00400     CHECK_FREAD(ret_val, "reading size of index array");
00401     CHECK_FEOF(ret_val, "reading size of index array");
00402     if (*reverseEndian) swap4_aligned(input_integer, rec_scale);
00403 
00404     /* when we have more then 2^30 atoms, tests like this one */
00405     /* are no longer meaningful and have to be bypassed...    */
00406     if (!hugefile && ((input_integer[0]+input_integer[1]) != ((*N)-(*NAMNF))*4L)) {
00407       return DCD_BADFORMAT;
00408     }
00409 
00410     ret_val = READ(fd, (*FREEINDEXES), ((ptrdiff_t) ((*N)-(*NAMNF)))*sizeof(int));
00411     CHECK_FREAD(ret_val, "reading size of index array");
00412     CHECK_FEOF(ret_val, "reading size of index array");
00413 
00414     if (*reverseEndian)
00415       swap4_aligned((*FREEINDEXES), ((*N)-(*NAMNF)));
00416 
00417     input_integer[1]=0;
00418     ret_val = READ(fd, input_integer, rec_scale*sizeof(int));
00419     CHECK_FREAD(ret_val, "reading size of index array");
00420     CHECK_FEOF(ret_val, "reading size of index array");
00421     if (*reverseEndian) swap4_aligned(input_integer, rec_scale);
00422 
00423     /* when we have more then 2^30 atoms, tests like this one */
00424     /* are no longer meaningful and have to be bypassed...    */
00425     if (!hugefile && ((input_integer[0]+input_integer[1]) != ((*N)-(*NAMNF))*4L)) {
00426       return DCD_BADFORMAT;
00427     }
00428   }
00429 
00430   return DCD_SUCCESS;
00431 }
00432 
00433 
00434 static int read_charmm_extrablock(fio_fd fd, int charmm, int reverseEndian,
00435                                   float *unitcell) {
00436   int i, input_integer[2], rec_scale;
00437 
00438   if (charmm & DCD_HAS_64BIT_REC) {
00439     rec_scale = RECSCALE64BIT;
00440   } else {
00441     rec_scale = RECSCALE32BIT;
00442   }
00443 
00444   if ((charmm & DCD_IS_CHARMM) && (charmm & DCD_HAS_EXTRA_BLOCK)) {
00445     /* Leading integer must be 48 */
00446     input_integer[1] = 0;
00447     if (fio_fread(input_integer, sizeof(int), rec_scale, fd) != rec_scale)
00448       return DCD_BADREAD; 
00449     if (reverseEndian) swap4_aligned(input_integer, rec_scale);
00450     if ((input_integer[0]+input_integer[1]) == 48) {
00451       double tmp[6];
00452       if (fio_fread(tmp, 48, 1, fd) != 1) return DCD_BADREAD;
00453       if (reverseEndian) 
00454         swap8_aligned(tmp, 6);
00455       for (i=0; i<6; i++) unitcell[i] = (float)tmp[i];
00456     } else {
00457       /* unrecognized block, just skip it */
00458       if (fio_fseek(fd, (input_integer[0]+input_integer[1]), FIO_SEEK_CUR)) return DCD_BADREAD;
00459     }
00460     if (fio_fread(input_integer, sizeof(int), rec_scale, fd) != rec_scale) return DCD_BADREAD; 
00461   } 
00462 
00463   return DCD_SUCCESS;
00464 }
00465 
00466 
00467 static int read_fixed_atoms(fio_fd fd, int N, int num_free, const int *indexes,
00468                             int reverseEndian, const float *fixedcoords, 
00469                             float *freeatoms, float *pos, int charmm) {
00470   int i, input_integer[2], rec_scale;
00471   
00472   if(charmm & DCD_HAS_64BIT_REC) {
00473     rec_scale=RECSCALE64BIT;
00474   } else {
00475     rec_scale=RECSCALE32BIT;
00476   }
00477   
00478   /* Read leading integer */
00479   input_integer[1]=0;
00480   if (fio_fread(input_integer, sizeof(int), rec_scale, fd) != rec_scale) return DCD_BADREAD;
00481   if (reverseEndian) swap4_aligned(input_integer, rec_scale);
00482   if ((input_integer[0]+input_integer[1]) != 4L*num_free) return DCD_BADFORMAT;
00483   
00484   /* Read free atom coordinates */
00485   if (fio_fread(freeatoms, 4L*num_free, 1, fd) != 1) return DCD_BADREAD;
00486   if (reverseEndian)
00487     swap4_aligned(freeatoms, num_free);
00488 
00489   /* Copy fixed and free atom coordinates into position buffer */
00490   memcpy(pos, fixedcoords, 4L*N);
00491   for (i=0; i<num_free; i++)
00492     pos[indexes[i]-1] = freeatoms[i];
00493 
00494   /* Read trailing integer */ 
00495   input_integer[1]=0;
00496   if (fio_fread(input_integer, sizeof(int), rec_scale, fd) != rec_scale) return DCD_BADREAD;
00497   if (reverseEndian) swap4_aligned(input_integer, rec_scale);
00498   if ((input_integer[0]+input_integer[1]) != 4L*num_free) return DCD_BADFORMAT;
00499 
00500   return DCD_SUCCESS;
00501 }
00502  
00503  
00504 static int read_charmm_4dim(fio_fd fd, int charmm, int reverseEndian) {
00505   int input_integer[2], rec_scale;
00506 
00507   if (charmm & DCD_HAS_64BIT_REC) {
00508     rec_scale=RECSCALE64BIT;
00509   } else {
00510     rec_scale=RECSCALE32BIT;
00511   }
00512     
00513   /* If this is a CHARMm file and contains a 4th dimension block, */
00514   /* we must skip past it to avoid problems                       */
00515   if ((charmm & DCD_IS_CHARMM) && (charmm & DCD_HAS_4DIMS)) {
00516     input_integer[1]=0;
00517     if (fio_fread(input_integer, sizeof(int), rec_scale, fd) != rec_scale) return DCD_BADREAD;  
00518     if (reverseEndian) swap4_aligned(input_integer, rec_scale);
00519     if (fio_fseek(fd, (input_integer[0]+input_integer[1]), FIO_SEEK_CUR)) return DCD_BADREAD;
00520     if (fio_fread(input_integer, sizeof(int), rec_scale, fd) != rec_scale) return DCD_BADREAD;  
00521   }
00522 
00523   return DCD_SUCCESS;
00524 }
00525 
00526 
00527 /* 
00528  * Read a dcd timestep from a dcd file
00529  * Input: fd - a file struct opened for binary reading, from which the 
00530  *             header information has already been read.
00531  *        natoms, nfixed, first, *freeind, reverse, charmm - the corresponding 
00532  *             items as set by read_dcdheader
00533  *        first - true if this is the first frame we are reading.
00534  *        x, y, z: space for natoms each of floats.
00535  *        unitcell - space for six floats to hold the unit cell data.  
00536  *                   Not set if no unit cell data is present.
00537  * Output: 0 on success, negative error code on failure.
00538  * Side effects: x, y, z contain the coordinates for the timestep read.
00539  *               unitcell holds unit cell data if present.
00540  */
00541 static int read_dcdstep(fio_fd fd, int N, float *X, float *Y, float *Z, 
00542                         float *unitcell, int num_fixed,
00543                         int first, int *indexes, float *fixedcoords, 
00544                         int reverseEndian, int charmm) {
00545   int ret_val;    /* Return value from read */
00546   ptrdiff_t rec_scale;
00547   int hugefile = (N > (1L<<30)) ? 1 : 0;
00548   int check_reclen = 1; /* Enable Fortran record length value safety checks */
00549 
00550   /* Fortran record length checks disabled for huge files or user request */
00551   if (hugefile || (getenv("VMDDCDNOCHECKRECLEN") != NULL))
00552     check_reclen = 0;
00553  
00554   if (charmm & DCD_HAS_64BIT_REC) {
00555     rec_scale=RECSCALE64BIT;
00556   } else {
00557     rec_scale=RECSCALE32BIT;
00558   }
00559   
00560   if ((num_fixed==0) || first) {
00561     /* temp storage for reading formatting info */
00562     /* note: has to be max size we'll ever use  */
00563     int tmpbuf[6L*RECSCALEMAX]; 
00564 
00565     fio_iovec iov[7];   /* I/O vector for fio_readv() call          */
00566     int i;
00567 
00568     /* if there are no fixed atoms or this is the first timestep read */
00569     /* then we read all coordinates normally.                         */
00570 
00571     /* read the charmm periodic cell information */
00572     /* XXX this too should be read together with the other items in a */
00573     /*     single fio_readv() call in order to prevent lots of extra  */
00574     /*     kernel/user context switches.                              */
00575     ret_val = read_charmm_extrablock(fd, charmm, reverseEndian, unitcell);
00576     if (ret_val) return ret_val;
00577 
00578     /* setup the I/O vector for the call to fio_readv() */
00579     iov[0].iov_base = (fio_caddr_t) &tmpbuf[0]; /* read format integer    */
00580     iov[0].iov_len  = rec_scale*sizeof(int);
00581 
00582     iov[1].iov_base = (fio_caddr_t) X;          /* read X coordinates     */
00583     iov[1].iov_len  = sizeof(float)*N;
00584 
00585     iov[2].iov_base = (fio_caddr_t) &tmpbuf[1*rec_scale]; /* read 2 format integers */
00586     iov[2].iov_len  = rec_scale*sizeof(int) * 2L;
00587 
00588     iov[3].iov_base = (fio_caddr_t) Y;          /* read Y coordinates     */
00589     iov[3].iov_len  = sizeof(float)*N;
00590 
00591     iov[4].iov_base = (fio_caddr_t) &tmpbuf[3L*rec_scale]; /* read 2 format integers */
00592     iov[4].iov_len  = rec_scale*sizeof(int) * 2L;
00593 
00594     iov[5].iov_base = (fio_caddr_t) Z;          /* read Y coordinates     */
00595     iov[5].iov_len  = sizeof(float)*N;
00596 
00597     iov[6].iov_base = (fio_caddr_t) &tmpbuf[5L*rec_scale]; /* read format integer    */
00598     iov[6].iov_len  = rec_scale*sizeof(int);
00599 
00600 #if 1
00601     /* Use fall-back code instead of readv():                            */
00602     /*  Some platforms implement readv() as user level code in libc,     */
00603     /*  and due to POSIX atomicity requirements for readv()/writev(),    */
00604     /*  they may copy data to internal temp buffers, which can kill      */
00605     /*  performance, and in cases when doing single I/O ops on large,    */
00606     /*  buffers, e.g. > 2GB, can fail with shorts reads or writes...     */
00607     /*  On such platforms it is best to avoid using readv()/writev()...  */
00608     {
00609       int readcnt = 0;
00610       readcnt =  fio_fread(iov[0].iov_base, iov[0].iov_len, 1, fd);
00611       readcnt += fio_fread(iov[1].iov_base, iov[1].iov_len, 1, fd);
00612       readcnt += fio_fread(iov[2].iov_base, iov[2].iov_len, 1, fd);
00613       readcnt += fio_fread(iov[3].iov_base, iov[3].iov_len, 1, fd);
00614       readcnt += fio_fread(iov[4].iov_base, iov[4].iov_len, 1, fd);
00615       readcnt += fio_fread(iov[5].iov_base, iov[5].iov_len, 1, fd);
00616       readcnt += fio_fread(iov[6].iov_base, iov[6].iov_len, 1, fd);
00617 
00618       /* if all records read correctly, then the reads are okay */
00619       if (readcnt != 7)
00620         return DCD_BADREAD;
00621     }
00622 #else
00623     /* check number of bytes actually read            */
00624     if (fio_readv(fd, &iov[0], 7) != ((fio_size_t) (rec_scale*6L*sizeof(int) + 3L*N*sizeof(float))))
00625       return DCD_BADREAD;
00626 #endif
00627 
00628     /* convert endianism if necessary */
00629     if (reverseEndian) {
00630       swap4_aligned(&tmpbuf[0], rec_scale*6L);
00631       swap4_aligned(X, N);
00632       swap4_aligned(Y, N);
00633       swap4_aligned(Z, N);
00634     }
00635 
00636     /* when we have more then 2^30 atoms, tests like this one */
00637     /* are no longer meaningful and have to be bypassed...    */
00638     if (check_reclen) {
00639       /* double-check the fortran format size values for safety */
00640       if (rec_scale == 1) {
00641         for (i=0; i<6; i++) {
00642           if (tmpbuf[i] != sizeof(float)*N) return DCD_BADFORMAT;
00643         }
00644       } else {
00645         for (i=0; i<6; i++) {
00646           if ((tmpbuf[2L*i]+tmpbuf[2L*i+1L]) != sizeof(float)*N) return DCD_BADFORMAT;
00647         }
00648       }
00649     }
00650 
00651     /* copy fixed atom coordinates into fixedcoords array if this was the */
00652     /* first timestep, to be used from now on.  We just copy all atoms.   */
00653     if (num_fixed && first) {
00654       memcpy(fixedcoords, X, N*sizeof(float));
00655       memcpy(fixedcoords+N, Y, N*sizeof(float));
00656       memcpy(fixedcoords+2L*N, Z, N*sizeof(float));
00657     }
00658 
00659     /* read in the optional charmm 4th array */
00660     /* XXX this too should be read together with the other items in a */
00661     /*     single fio_readv() call in order to prevent lots of extra  */
00662     /*     kernel/user context switches.                              */
00663     ret_val = read_charmm_4dim(fd, charmm, reverseEndian);
00664     if (ret_val) return ret_val;
00665   } else {
00666     /* if there are fixed atoms, and this isn't the first frame, then we */
00667     /* only read in the non-fixed atoms for all subsequent timesteps.    */
00668     ret_val = read_charmm_extrablock(fd, charmm, reverseEndian, unitcell);
00669     if (ret_val) return ret_val;
00670     ret_val = read_fixed_atoms(fd, N, N-num_fixed, indexes, reverseEndian,
00671                                fixedcoords, fixedcoords+3L*N, X, charmm);
00672     if (ret_val) return ret_val;
00673     ret_val = read_fixed_atoms(fd, N, N-num_fixed, indexes, reverseEndian,
00674                                fixedcoords+N, fixedcoords+3L*N, Y, charmm);
00675     if (ret_val) return ret_val;
00676     ret_val = read_fixed_atoms(fd, N, N-num_fixed, indexes, reverseEndian,
00677                                fixedcoords+2*N, fixedcoords+3L*N, Z, charmm);
00678     if (ret_val) return ret_val;
00679     ret_val = read_charmm_4dim(fd, charmm, reverseEndian);
00680     if (ret_val) return ret_val;
00681   }
00682 
00683   return DCD_SUCCESS;
00684 }
00685 
00686 
00687 /* 
00688  * Skip past a timestep.  If there are fixed atoms, this cannot be used with
00689  * the first timestep.  
00690  * Input: fd - a file struct from which the header has already been read
00691  *        natoms - number of atoms per timestep
00692  *        nfixed - number of fixed atoms
00693  *        charmm - charmm flags as returned by read_dcdheader
00694  * Output: 0 on success, negative error code on failure.
00695  * Side effects: One timestep will be skipped; fd will be positioned at the
00696  *               next timestep.
00697  */
00698 static int skip_dcdstep(fio_fd fd, int natoms, int nfixed, int charmm) {
00699   ptrdiff_t seekoffset = 0;
00700   ptrdiff_t rec_scale;
00701 
00702   if (charmm & DCD_HAS_64BIT_REC) {
00703     rec_scale=RECSCALE64BIT;
00704   } else {
00705     rec_scale=RECSCALE32BIT;
00706   }
00707 
00708   /* Skip charmm extra block */
00709   if ((charmm & DCD_IS_CHARMM) && (charmm & DCD_HAS_EXTRA_BLOCK)) {
00710     seekoffset += 4L*rec_scale + 48L + 4L*rec_scale;
00711   }
00712 
00713   /* For each atom set, seek past an int, the free atoms, and another int. */
00714   seekoffset += 3L * (2L*rec_scale + natoms - nfixed) * 4L;
00715 
00716   /* Assume that charmm 4th dim is the same size as the other three. */
00717   if ((charmm & DCD_IS_CHARMM) && (charmm & DCD_HAS_4DIMS)) {
00718     seekoffset += (2L*rec_scale + natoms - nfixed) * 4L;
00719   }
00720  
00721   if (fio_fseek(fd, seekoffset, FIO_SEEK_CUR)) return DCD_BADEOF;
00722 
00723   return DCD_SUCCESS;
00724 }
00725 
00726 
00727 /* 
00728  * Write a timestep to a dcd file
00729  * Input: fd - a file struct for which a dcd header has already been written
00730  *       curframe: Count of frames written to this file, starting with 1.
00731  *        curstep: Count of timesteps elapsed = istart + curframe * nsavc.
00732  *         natoms: number of elements in x, y, z arrays
00733  *        x, y, z: pointers to atom coordinates
00734  * Output: 0 on success, negative error code on failure.
00735  * Side effects: coordinates are written to the dcd file.
00736  */
00737 static int write_dcdstep(fio_fd fd, int curframe, int curstep, int N, 
00738                   const float *X, const float *Y, const float *Z, 
00739                   const double *unitcell, int charmm) {
00740   int out_integer;
00741 
00742   if (charmm) {
00743     /* write out optional unit cell */
00744     if (unitcell != NULL) {
00745       out_integer = 48; /* 48 bytes (6 floats) */
00746       fio_write_int32(fd, out_integer);
00747       WRITE(fd, unitcell, out_integer);
00748       fio_write_int32(fd, out_integer);
00749     }
00750   }
00751 
00752   /* write out coordinates */
00753   out_integer = N*4; /* N*4 bytes per X/Y/Z array (N floats per array) */
00754   fio_write_int32(fd, out_integer);
00755   if (fio_fwrite((void *) X, out_integer, 1, fd) != 1) return DCD_BADWRITE;
00756   fio_write_int32(fd, out_integer);
00757   fio_write_int32(fd, out_integer);
00758   if (fio_fwrite((void *) Y, out_integer, 1, fd) != 1) return DCD_BADWRITE;
00759   fio_write_int32(fd, out_integer);
00760   fio_write_int32(fd, out_integer);
00761   if (fio_fwrite((void *) Z, out_integer, 1, fd) != 1) return DCD_BADWRITE;
00762   fio_write_int32(fd, out_integer);
00763 
00764   /* update the DCD header information */
00765   fio_fseek(fd, NFILE_POS, FIO_SEEK_SET);
00766   fio_write_int32(fd, curframe);
00767   fio_fseek(fd, NSTEP_POS, FIO_SEEK_SET);
00768   fio_write_int32(fd, curstep);
00769   fio_fseek(fd, 0, FIO_SEEK_END);
00770 
00771   return DCD_SUCCESS;
00772 }
00773 
00774 
00775 /*
00776  * Write a header for a new dcd file
00777  * Input: fd - file struct opened for binary writing
00778  *        remarks - string to be put in the remarks section of the header.  
00779  *                  The string will be truncated to 70 characters.
00780  *        natoms, istart, nsavc, delta - see comments in read_dcdheader
00781  * Output: 0 on success, negative error code on failure.
00782  * Side effects: Header information is written to the dcd file.
00783  */
00784 static int write_dcdheader(fio_fd fd, const char *remarks, int N, 
00785                     int ISTART, int NSAVC, double DELTA, int with_unitcell,
00786                     int charmm) {
00787   int out_integer;
00788   float out_float;
00789   char title_string[200];
00790   time_t cur_time;
00791   struct tm *tmbuf;
00792   char time_str[81];
00793 
00794   out_integer = 84;
00795   WRITE(fd, (char *) & out_integer, sizeof(int));
00796   strcpy(title_string, "CORD");
00797   WRITE(fd, title_string, 4);
00798   fio_write_int32(fd, 0);      /* Number of frames in file, none written yet   */
00799   fio_write_int32(fd, ISTART); /* Starting timestep                            */
00800   fio_write_int32(fd, NSAVC);  /* Timesteps between frames written to the file */
00801   fio_write_int32(fd, 0);      /* Number of timesteps in simulation            */
00802   fio_write_int32(fd, 0);      /* NAMD writes NSTEP or ISTART - NSAVC here?    */
00803   fio_write_int32(fd, 0);
00804   fio_write_int32(fd, 0);
00805   fio_write_int32(fd, 0);
00806   fio_write_int32(fd, 0);
00807   if (charmm) {
00808     out_float = DELTA;
00809     WRITE(fd, (char *) &out_float, sizeof(float));
00810     if (with_unitcell) {
00811       fio_write_int32(fd, 1);
00812     } else {
00813       fio_write_int32(fd, 0);
00814     }
00815   } else {
00816     WRITE(fd, (char *) &DELTA, sizeof(double));
00817   }
00818   fio_write_int32(fd, 0);
00819   fio_write_int32(fd, 0);
00820   fio_write_int32(fd, 0);
00821   fio_write_int32(fd, 0);
00822   fio_write_int32(fd, 0);
00823   fio_write_int32(fd, 0);
00824   fio_write_int32(fd, 0);
00825   fio_write_int32(fd, 0);
00826   if (charmm) {
00827     fio_write_int32(fd, 24); /* Pretend to be CHARMM version 24 */
00828   } else {
00829     fio_write_int32(fd, 0);
00830   }
00831   fio_write_int32(fd, 84);
00832   fio_write_int32(fd, 164);
00833   fio_write_int32(fd, 2);
00834 
00835   strncpy(title_string, remarks, 80);
00836   title_string[79] = '\0';
00837   WRITE(fd, title_string, 80);
00838 
00839   cur_time=time(NULL);
00840   tmbuf=localtime(&cur_time);
00841   strftime(time_str, 80, "REMARKS Created %d %B, %Y at %R", tmbuf);
00842   WRITE(fd, time_str, 80);
00843 
00844   fio_write_int32(fd, 164);
00845   fio_write_int32(fd, 4);
00846   fio_write_int32(fd, N);
00847   fio_write_int32(fd, 4);
00848 
00849   return DCD_SUCCESS;
00850 }
00851 
00852 
00853 /*
00854  * clean up dcd data
00855  * Input: nfixed, freeind - elements as returned by read_dcdheader
00856  * Output: None
00857  * Side effects: Space pointed to by freeind is freed if necessary.
00858  */
00859 static void close_dcd_read(int *indexes, float *fixedcoords) {
00860   free(indexes);
00861   free(fixedcoords);
00862 }
00863 
00864 
00865 static void *open_dcd_read(const char *path, const char *filetype, 
00866     int *natoms) {
00867   dcdhandle *dcd;
00868   fio_fd fd;
00869   int rc;
00870   struct stat stbuf;
00871 
00872   if (!path) return NULL;
00873 
00874 #if !(defined(_MSC_VER) && defined(FASTIO_NATIVEWIN32))
00875   /* See if the file exists, and get its size */
00876   memset(&stbuf, 0, sizeof(struct stat));
00877   if (stat(path, &stbuf)) {
00878     printf("dcdplugin) Could not access file '%s'.\n", path);
00879     return NULL;
00880   }
00881 #endif
00882 
00883   if (fio_open(path, FIO_READ, &fd) < 0) {
00884     printf("dcdplugin) Could not open file '%s' for reading.\n", path);
00885     return NULL;
00886   }
00887 
00888   dcd = (dcdhandle *)malloc(sizeof(dcdhandle));
00889   memset(dcd, 0, sizeof(dcdhandle));
00890   dcd->fd = fd;
00891 
00892   if ((rc = read_dcdheader(dcd->fd, &dcd->natoms, &dcd->nsets, &dcd->istart, 
00893          &dcd->nsavc, &dcd->delta, &dcd->nfixed, &dcd->freeind, 
00894          &dcd->fixedcoords, &dcd->reverse, &dcd->charmm))) {
00895     print_dcderror("read_dcdheader", rc);
00896     fio_fclose(dcd->fd);
00897     free(dcd);
00898     return NULL;
00899   }
00900 
00901   /*
00902    * Check that the file is big enough to really hold the number of sets
00903    * it claims to have.  Then we'll use nsets to keep track of where EOF
00904    * should be.
00905    */
00906   {
00907     fio_size_t ndims, firstframesize, framesize, extrablocksize;
00908     fio_size_t trjsize, filesize, curpos;
00909     int newnsets;
00910 
00911     extrablocksize = dcd->charmm & DCD_HAS_EXTRA_BLOCK ? 48 + 8 : 0;
00912     ndims = dcd->charmm & DCD_HAS_4DIMS ? 4 : 3;
00913     firstframesize = (dcd->natoms+2) * ndims * sizeof(float) + extrablocksize;
00914     framesize = (dcd->natoms-dcd->nfixed+2) * ndims * sizeof(float) 
00915       + extrablocksize;
00916 
00917     /* 
00918      * It's safe to use ftell, even though ftell returns a long, because the 
00919      * header size is < 4GB.
00920      */
00921 
00922     curpos = fio_ftell(dcd->fd); /* save current offset (end of header) */
00923 
00924 #if defined(_MSC_VER) && defined(FASTIO_NATIVEWIN32)
00925     /* the stat() call is not 64-bit savvy on Windows             */
00926     /* so we have to use the fastio fseek/ftell routines for this */
00927     /* until we add a portable filesize routine for this purpose  */
00928     fio_fseek(dcd->fd, 0, FIO_SEEK_END);       /* seek to end of file */
00929     filesize = fio_ftell(dcd->fd);
00930     fio_fseek(dcd->fd, curpos, FIO_SEEK_SET);  /* return to end of header */
00931 #else
00932     filesize = stbuf.st_size; /* this works ok on Unix machines */
00933 #endif
00934     trjsize = filesize - curpos - firstframesize;
00935     if (trjsize < 0) {
00936       printf("dcdplugin) file '%s' appears to contain no timesteps.\n", path);
00937       fio_fclose(dcd->fd);
00938       free(dcd);
00939       return NULL;
00940     }
00941 
00942     newnsets = trjsize / framesize + 1;
00943 
00944     if (dcd->nsets > 0 && newnsets != dcd->nsets) {
00945       printf("dcdplugin) Warning: DCD header claims %d frames, but \n"
00946              "dcdplugin) file size (%ld) indicates there are actually \n"
00947              "%d frames of size (%ld)\n", 
00948              dcd->nsets, trjsize, newnsets, framesize);
00949     }
00950 
00951     dcd->nsets = newnsets; 
00952     dcd->setsread = 0;
00953   }
00954 
00955   dcd->first = 1;
00956   dcd->x = (float *)malloc(dcd->natoms * sizeof(float));
00957   dcd->y = (float *)malloc(dcd->natoms * sizeof(float));
00958   dcd->z = (float *)malloc(dcd->natoms * sizeof(float));
00959   if (!dcd->x || !dcd->y || !dcd->z) {
00960     printf("dcdplugin) Unable to allocate space for %d atoms.\n", dcd->natoms);
00961     if (dcd->x)
00962       free(dcd->x);
00963     if (dcd->y)
00964       free(dcd->y);
00965     if (dcd->z)
00966       free(dcd->z);
00967     fio_fclose(dcd->fd);
00968     free(dcd);
00969     return NULL;
00970   }
00971   *natoms = dcd->natoms;
00972   return dcd;
00973 }
00974 
00975 
00976 static int read_next_timestep(void *v, int natoms, molfile_timestep_t *ts) {
00977   dcdhandle *dcd;
00978   int i, j, rc;
00979   float unitcell[6];
00980   unitcell[0] = unitcell[2] = unitcell[5] = 0.0f;
00981   unitcell[1] = unitcell[3] = unitcell[4] = 90.0f;
00982   dcd = (dcdhandle *)v;
00983 
00984   /* Check for EOF here; that way all EOF's encountered later must be errors */
00985   if (dcd->setsread == dcd->nsets) return MOLFILE_EOF;
00986   dcd->setsread++;
00987   if (!ts) {
00988     if (dcd->first && dcd->nfixed) {
00989       /* We can't just skip it because we need the fixed atom coordinates */
00990       rc = read_dcdstep(dcd->fd, dcd->natoms, dcd->x, dcd->y, dcd->z, 
00991           unitcell, dcd->nfixed, dcd->first, dcd->freeind, dcd->fixedcoords, 
00992              dcd->reverse, dcd->charmm);
00993       dcd->first = 0;
00994       return rc; /* XXX this needs to be updated */
00995     }
00996     dcd->first = 0;
00997     /* XXX this needs to be changed */
00998     return skip_dcdstep(dcd->fd, dcd->natoms, dcd->nfixed, dcd->charmm);
00999   }
01000   rc = read_dcdstep(dcd->fd, dcd->natoms, dcd->x, dcd->y, dcd->z, unitcell,
01001              dcd->nfixed, dcd->first, dcd->freeind, dcd->fixedcoords, 
01002              dcd->reverse, dcd->charmm);
01003   dcd->first = 0;
01004   if (rc < 0) {  
01005     print_dcderror("read_dcdstep", rc);
01006     return MOLFILE_ERROR;
01007   }
01008 
01009   /* copy timestep data from plugin-local buffers to VMD's buffer */
01010   /* XXX 
01011    *   This code is still the root of all evil.  Just doing this extra copy
01012    *   cuts the I/O rate of the DCD reader from 728 MB/sec down to
01013    *   394 MB/sec when reading from a ram filesystem.  
01014    *   For a physical disk filesystem, the I/O rate goes from 
01015    *   187 MB/sec down to 122 MB/sec.  Clearly this extra copy has to go.
01016    */
01017   {
01018     int natoms = dcd->natoms;
01019     float *nts = ts->coords;
01020     const float *bufx = dcd->x;
01021     const float *bufy = dcd->y;
01022     const float *bufz = dcd->z;
01023 
01024     for (i=0, j=0; i<natoms; i++, j+=3) {
01025       nts[j    ] = bufx[i];
01026       nts[j + 1] = bufy[i];
01027       nts[j + 2] = bufz[i];
01028     }
01029   }
01030 
01031   ts->A = unitcell[0];
01032   ts->B = unitcell[2];
01033   ts->C = unitcell[5];
01034 
01035   if (unitcell[1] >= -1.0 && unitcell[1] <= 1.0 &&
01036       unitcell[3] >= -1.0 && unitcell[3] <= 1.0 &&
01037       unitcell[4] >= -1.0 && unitcell[4] <= 1.0) {
01038     /* This file was generated by CHARMM, or by NAMD > 2.5, with the angle */
01039     /* cosines of the periodic cell angles written to the DCD file.        */ 
01040     /* This formulation improves rounding behavior for orthogonal cells    */
01041     /* so that the angles end up at precisely 90 degrees, unlike acos().   */
01042     ts->alpha = 90.0 - asin(unitcell[4]) * 90.0 / M_PI_2; /* cosBC */
01043     ts->beta  = 90.0 - asin(unitcell[3]) * 90.0 / M_PI_2; /* cosAC */
01044     ts->gamma = 90.0 - asin(unitcell[1]) * 90.0 / M_PI_2; /* cosAB */
01045   } else {
01046     /* This file was likely generated by NAMD 2.5 and the periodic cell    */
01047     /* angles are specified in degrees rather than angle cosines.          */
01048     ts->alpha = unitcell[4]; /* angle between B and C */
01049     ts->beta  = unitcell[3]; /* angle between A and C */
01050     ts->gamma = unitcell[1]; /* angle between A and B */
01051   }
01052  
01053   return MOLFILE_SUCCESS;
01054 }
01055  
01056 
01057 static void close_file_read(void *v) {
01058   dcdhandle *dcd = (dcdhandle *)v;
01059   close_dcd_read(dcd->freeind, dcd->fixedcoords);
01060   fio_fclose(dcd->fd);
01061   free(dcd->x);
01062   free(dcd->y);
01063   free(dcd->z);
01064   free(dcd); 
01065 }
01066 
01067 
01068 static void *open_dcd_write(const char *path, const char *filetype, 
01069     int natoms) {
01070   dcdhandle *dcd;
01071   fio_fd fd;
01072   int rc;
01073   int istart, nsavc;
01074   double delta;
01075   int with_unitcell;
01076   int charmm;
01077 
01078   if (fio_open(path, FIO_WRITE, &fd) < 0) {
01079     printf("dcdplugin) Could not open file '%s' for writing\n", path);
01080     return NULL;
01081   }
01082 
01083   dcd = (dcdhandle *)malloc(sizeof(dcdhandle));
01084   memset(dcd, 0, sizeof(dcdhandle));
01085   dcd->fd = fd;
01086 
01087   istart = 0;             /* starting timestep of DCD file                  */
01088   nsavc = 1;              /* number of timesteps between written DCD frames */
01089   delta = 1.0;            /* length of a timestep                           */
01090 
01091   if (getenv("VMDDCDWRITEXPLORFORMAT") != NULL) {
01092     with_unitcell = 0;      /* no unit cell info */
01093     charmm = DCD_IS_XPLOR;  /* X-PLOR format */
01094     printf("dcdplugin) WARNING: Writing DCD file in X-PLOR format, \n");
01095     printf("dcdplugin) WARNING: unit cell information will be lost!\n");
01096   } else {
01097     with_unitcell = 1;      /* contains unit cell infor (Charmm format) */
01098     charmm = DCD_IS_CHARMM; /* charmm-formatted DCD file                */ 
01099     if (with_unitcell) 
01100       charmm |= DCD_HAS_EXTRA_BLOCK;
01101   }
01102  
01103   rc = write_dcdheader(dcd->fd, "Created by DCD plugin", natoms, 
01104                        istart, nsavc, delta, with_unitcell, charmm);
01105 
01106   if (rc < 0) {
01107     print_dcderror("write_dcdheader", rc);
01108     fio_fclose(dcd->fd);
01109     free(dcd);
01110     return NULL;
01111   }
01112 
01113   dcd->natoms = natoms;
01114   dcd->nsets = 0;
01115   dcd->istart = istart;
01116   dcd->nsavc = nsavc;
01117   dcd->with_unitcell = with_unitcell;
01118   dcd->charmm = charmm;
01119   dcd->x = (float *)malloc(natoms * sizeof(float));
01120   dcd->y = (float *)malloc(natoms * sizeof(float));
01121   dcd->z = (float *)malloc(natoms * sizeof(float));
01122   return dcd;
01123 }
01124 
01125 
01126 static int write_timestep(void *v, const molfile_timestep_t *ts) { 
01127   dcdhandle *dcd = (dcdhandle *)v;
01128   int i, rc, curstep;
01129   float *pos = ts->coords;
01130   double unitcell[6];
01131   unitcell[0] = unitcell[2] = unitcell[5] = 0.0f;
01132   unitcell[1] = unitcell[3] = unitcell[4] = 90.0f;
01133 
01134   /* copy atom coords into separate X/Y/Z arrays for writing */
01135   for (i=0; i<dcd->natoms; i++) {
01136     dcd->x[i] = *(pos++); 
01137     dcd->y[i] = *(pos++); 
01138     dcd->z[i] = *(pos++); 
01139   }
01140   dcd->nsets++;
01141   curstep = dcd->istart + dcd->nsets * dcd->nsavc;
01142 
01143   unitcell[0] = ts->A;
01144   unitcell[2] = ts->B;
01145   unitcell[5] = ts->C;
01146   unitcell[1] = sin((M_PI_2 / 90.0) * (90.0 - ts->gamma)); /* cosAB */
01147   unitcell[3] = sin((M_PI_2 / 90.0) * (90.0 - ts->beta));  /* cosAC */
01148   unitcell[4] = sin((M_PI_2 / 90.0) * (90.0 - ts->alpha)); /* cosBC */
01149 
01150   rc = write_dcdstep(dcd->fd, dcd->nsets, curstep, dcd->natoms, 
01151                      dcd->x, dcd->y, dcd->z,
01152                      dcd->with_unitcell ? unitcell : NULL,
01153                      dcd->charmm);
01154   if (rc < 0) {
01155     print_dcderror("write_dcdstep", rc);
01156     return MOLFILE_ERROR;
01157   }
01158 
01159   return MOLFILE_SUCCESS;
01160 }
01161 
01162 static void close_file_write(void *v) {
01163   dcdhandle *dcd = (dcdhandle *)v;
01164   fio_fclose(dcd->fd);
01165   free(dcd->x);
01166   free(dcd->y);
01167   free(dcd->z);
01168   free(dcd);
01169 }
01170 
01171 
01172 /*
01173  * Initialization stuff here
01174  */
01175 static molfile_plugin_t plugin;
01176 
01177 VMDPLUGIN_API int VMDPLUGIN_init() {
01178   memset(&plugin, 0, sizeof(molfile_plugin_t));
01179   plugin.abiversion = vmdplugin_ABIVERSION;
01180   plugin.type = MOLFILE_PLUGIN_TYPE;
01181   plugin.name = "dcd";
01182   plugin.prettyname = "CHARMM,NAMD,XPLOR DCD Trajectory";
01183   plugin.author = "Axel Kohlmeyer, Justin Gullingsrud, John Stone";
01184   plugin.majorv = 1;
01185   plugin.minorv = 18;
01186   plugin.is_reentrant = VMDPLUGIN_THREADSAFE;
01187   plugin.filename_extension = "dcd";
01188   plugin.open_file_read = open_dcd_read;
01189   plugin.read_next_timestep = read_next_timestep;
01190   plugin.close_file_read = close_file_read;
01191   plugin.open_file_write = open_dcd_write;
01192   plugin.write_timestep = write_timestep;
01193   plugin.close_file_write = close_file_write;
01194   return VMDPLUGIN_SUCCESS;
01195 }
01196 
01197 VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
01198   (*cb)(v, (vmdplugin_t *)&plugin);
01199   return VMDPLUGIN_SUCCESS;
01200 }
01201 
01202 VMDPLUGIN_API int VMDPLUGIN_fini() {
01203   return VMDPLUGIN_SUCCESS;
01204 }
01205 
01206   
01207 #ifdef TEST_DCDPLUGIN
01208 
01209 #include <sys/time.h>
01210 
01211 /* get the time of day from the system clock, and store it (in seconds) */
01212 double time_of_day(void) {
01213 #if defined(_MSC_VER)
01214   double t;
01215 
01216   t = GetTickCount();
01217   t = t / 1000.0;
01218 
01219   return t;
01220 #else
01221   struct timeval tm;
01222   struct timezone tz;
01223 
01224   gettimeofday(&tm, &tz);
01225   return((double)(tm.tv_sec) + (double)(tm.tv_usec)/1000000.0);
01226 #endif
01227 }
01228 
01229 int main(int argc, char *argv[]) {
01230   molfile_timestep_t timestep;
01231   void *v;
01232   dcdhandle *dcd;
01233   int i, natoms;
01234   float sizeMB =0.0, totalMB = 0.0;
01235   double starttime, endtime, totaltime = 0.0;
01236 
01237   while (--argc) {
01238     ++argv; 
01239     natoms = 0;
01240     v = open_dcd_read(*argv, "dcd", &natoms);
01241     if (!v) {
01242       fprintf(stderr, "main) open_dcd_read failed for file %s\n", *argv);
01243       return 1;
01244     }
01245     dcd = (dcdhandle *)v;
01246     sizeMB = ((natoms * 3.0) * dcd->nsets * 4.0) / (1024.0 * 1024.0);
01247     totalMB += sizeMB; 
01248     printf("main) file: %s\n", *argv);
01249     printf("  %d atoms, %d frames, size: %6.1fMB\n", natoms, dcd->nsets, sizeMB);
01250 
01251     starttime = time_of_day();
01252     timestep.coords = (float *)malloc(3*sizeof(float)*natoms);
01253     for (i=0; i<dcd->nsets; i++) {
01254       int rc = read_next_timestep(v, natoms, &timestep);
01255       if (rc) {
01256         fprintf(stderr, "error in read_next_timestep on frame %d\n", i);
01257         return 1;
01258       }
01259     }
01260     endtime = time_of_day();
01261     close_file_read(v);
01262     totaltime += endtime - starttime;
01263     printf("  Time: %5.1f seconds\n", endtime - starttime);
01264     printf("  Speed: %5.1f MB/sec, %5.1f timesteps/sec\n", sizeMB / (endtime - starttime), (dcd->nsets / (endtime - starttime)));
01265   }
01266   printf("Overall Size: %6.1f MB\n", totalMB);
01267   printf("Overall Time: %6.1f seconds\n", totaltime);
01268   printf("Overall Speed: %5.1f MB/sec\n", totalMB / totaltime);
01269   return 0;
01270 }
01271       
01272 #endif
01273 

Generated on Thu Nov 7 03:09:35 2024 for VMD Plugins (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002