Main Page   Alphabetical List   Compound List   File List   Compound Members   File Members   Related Pages  

dcdplugin.c

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

Generated on Sat May 18 02:03:37 2013 for VMD Plugins (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002