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

Generated on Fri Jan 24 02:56:16 2020 for VMD Plugins (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002