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.77 $       $Date: 2013/06/04 21:12:12 $
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 #if 1
00577     /* Use fall-back code instead of readv():
00578     /*  Some platforms implement readv() as user level code in libc,     */
00579     /*  and due to POSIX atomicity requirements for readv()/writev(),    */
00580     /*  they may copy data to internal temp buffers, which can kill      */
00581     /*  performance, and in cases when doing single I/O ops on large,    */
00582     /*  buffers, e.g. > 2GB, can fail with shorts reads or writes...     */
00583     /*  On such platforms it is best to avoid using readv()/writev()...  */
00584     {
00585       int readcnt = 0;
00586       readlen = 0;
00587       readcnt =  fio_fread(iov[0].iov_base, iov[0].iov_len, 1, fd);
00588       readcnt += fio_fread(iov[1].iov_base, iov[1].iov_len, 1, fd);
00589       readcnt += fio_fread(iov[2].iov_base, iov[2].iov_len, 1, fd);
00590       readcnt += fio_fread(iov[3].iov_base, iov[3].iov_len, 1, fd);
00591       readcnt += fio_fread(iov[4].iov_base, iov[4].iov_len, 1, fd);
00592       readcnt += fio_fread(iov[5].iov_base, iov[5].iov_len, 1, fd);
00593       readcnt += fio_fread(iov[6].iov_base, iov[6].iov_len, 1, fd);
00594 
00595       /* if both records read correctly, then the reads are okay */
00596       if (readcnt != 7)
00597         return DCD_BADREAD;
00598     }
00599 #else
00600     readlen = fio_readv(fd, &iov[0], 7);
00601     if (readlen != (rec_scale*6*sizeof(int) + 3*N*sizeof(float)))
00602       return DCD_BADREAD;
00603 #endif
00604 
00605     /* convert endianism if necessary */
00606     if (reverseEndian) {
00607       swap4_aligned(&tmpbuf[0], rec_scale*6);
00608       swap4_aligned(X, N);
00609       swap4_aligned(Y, N);
00610       swap4_aligned(Z, N);
00611     }
00612 
00613     /* double-check the fortran format size values for safety */
00614     if(rec_scale == 1) {
00615       for (i=0; i<6; i++) {
00616         if (tmpbuf[i] != sizeof(float)*N) return DCD_BADFORMAT;
00617       }
00618     } else {
00619       for (i=0; i<6; i++) {
00620           if ((tmpbuf[2*i]+tmpbuf[2*i+1]) != sizeof(float)*N) return DCD_BADFORMAT;
00621       }
00622     }
00623 
00624     /* copy fixed atom coordinates into fixedcoords array if this was the */
00625     /* first timestep, to be used from now on.  We just copy all atoms.   */
00626     if (num_fixed && first) {
00627       memcpy(fixedcoords, X, N*sizeof(float));
00628       memcpy(fixedcoords+N, Y, N*sizeof(float));
00629       memcpy(fixedcoords+2*N, Z, N*sizeof(float));
00630     }
00631 
00632     /* read in the optional charmm 4th array */
00633     /* XXX this too should be read together with the other items in a */
00634     /*     single fio_readv() call in order to prevent lots of extra  */
00635     /*     kernel/user context switches.                              */
00636     ret_val = read_charmm_4dim(fd, charmm, reverseEndian);
00637     if (ret_val) return ret_val;
00638   } else {
00639     /* if there are fixed atoms, and this isn't the first frame, then we */
00640     /* only read in the non-fixed atoms for all subsequent timesteps.    */
00641     ret_val = read_charmm_extrablock(fd, charmm, reverseEndian, unitcell);
00642     if (ret_val) return ret_val;
00643     ret_val = read_fixed_atoms(fd, N, N-num_fixed, indexes, reverseEndian,
00644                                fixedcoords, fixedcoords+3*N, X, charmm);
00645     if (ret_val) return ret_val;
00646     ret_val = read_fixed_atoms(fd, N, N-num_fixed, indexes, reverseEndian,
00647                                fixedcoords+N, fixedcoords+3*N, Y, charmm);
00648     if (ret_val) return ret_val;
00649     ret_val = read_fixed_atoms(fd, N, N-num_fixed, indexes, reverseEndian,
00650                                fixedcoords+2*N, fixedcoords+3*N, Z, charmm);
00651     if (ret_val) return ret_val;
00652     ret_val = read_charmm_4dim(fd, charmm, reverseEndian);
00653     if (ret_val) return ret_val;
00654   }
00655 
00656   return DCD_SUCCESS;
00657 }
00658 
00659 
00660 /* 
00661  * Skip past a timestep.  If there are fixed atoms, this cannot be used with
00662  * the first timestep.  
00663  * Input: fd - a file struct from which the header has already been read
00664  *        natoms - number of atoms per timestep
00665  *        nfixed - number of fixed atoms
00666  *        charmm - charmm flags as returned by read_dcdheader
00667  * Output: 0 on success, negative error code on failure.
00668  * Side effects: One timestep will be skipped; fd will be positioned at the
00669  *               next timestep.
00670  */
00671 static int skip_dcdstep(fio_fd fd, int natoms, int nfixed, int charmm) {
00672   
00673   int seekoffset = 0;
00674   int rec_scale;
00675 
00676   if (charmm & DCD_HAS_64BIT_REC) {
00677     rec_scale=RECSCALE64BIT;
00678   } else {
00679     rec_scale=RECSCALE32BIT;
00680   }
00681 
00682   /* Skip charmm extra block */
00683   if ((charmm & DCD_IS_CHARMM) && (charmm & DCD_HAS_EXTRA_BLOCK)) {
00684     seekoffset += 4*rec_scale + 48 + 4*rec_scale;
00685   }
00686 
00687   /* For each atom set, seek past an int, the free atoms, and another int. */
00688   seekoffset += 3 * (2*rec_scale + natoms - nfixed) * 4;
00689 
00690   /* Assume that charmm 4th dim is the same size as the other three. */
00691   if ((charmm & DCD_IS_CHARMM) && (charmm & DCD_HAS_4DIMS)) {
00692     seekoffset += (2*rec_scale + natoms - nfixed) * 4;
00693   }
00694  
00695   if (fio_fseek(fd, seekoffset, FIO_SEEK_CUR)) return DCD_BADEOF;
00696 
00697   return DCD_SUCCESS;
00698 }
00699 
00700 
00701 /* 
00702  * Write a timestep to a dcd file
00703  * Input: fd - a file struct for which a dcd header has already been written
00704  *       curframe: Count of frames written to this file, starting with 1.
00705  *       curstep: Count of timesteps elapsed = istart + curframe * nsavc.
00706  *        natoms - number of elements in x, y, z arrays
00707  *        x, y, z: pointers to atom coordinates
00708  * Output: 0 on success, negative error code on failure.
00709  * Side effects: coordinates are written to the dcd file.
00710  */
00711 static int write_dcdstep(fio_fd fd, int curframe, int curstep, int N, 
00712                   const float *X, const float *Y, const float *Z, 
00713                   const double *unitcell, int charmm) {
00714   int out_integer;
00715 
00716   if (charmm) {
00717     /* write out optional unit cell */
00718     if (unitcell != NULL) {
00719       out_integer = 48; /* 48 bytes (6 floats) */
00720       fio_write_int32(fd, out_integer);
00721       WRITE(fd, unitcell, out_integer);
00722       fio_write_int32(fd, out_integer);
00723     }
00724   }
00725 
00726   /* write out coordinates */
00727   out_integer = N*4; /* N*4 bytes per X/Y/Z array (N floats per array) */
00728   fio_write_int32(fd, out_integer);
00729   if (fio_fwrite((void *) X, out_integer, 1, fd) != 1) return DCD_BADWRITE;
00730   fio_write_int32(fd, out_integer);
00731   fio_write_int32(fd, out_integer);
00732   if (fio_fwrite((void *) Y, out_integer, 1, fd) != 1) return DCD_BADWRITE;
00733   fio_write_int32(fd, out_integer);
00734   fio_write_int32(fd, out_integer);
00735   if (fio_fwrite((void *) Z, out_integer, 1, fd) != 1) return DCD_BADWRITE;
00736   fio_write_int32(fd, out_integer);
00737 
00738   /* update the DCD header information */
00739   fio_fseek(fd, NFILE_POS, FIO_SEEK_SET);
00740   fio_write_int32(fd, curframe);
00741   fio_fseek(fd, NSTEP_POS, FIO_SEEK_SET);
00742   fio_write_int32(fd, curstep);
00743   fio_fseek(fd, 0, FIO_SEEK_END);
00744 
00745   return DCD_SUCCESS;
00746 }
00747 
00748 /*
00749  * Write a header for a new dcd file
00750  * Input: fd - file struct opened for binary writing
00751  *        remarks - string to be put in the remarks section of the header.  
00752  *                  The string will be truncated to 70 characters.
00753  *        natoms, istart, nsavc, delta - see comments in read_dcdheader
00754  * Output: 0 on success, negative error code on failure.
00755  * Side effects: Header information is written to the dcd file.
00756  */
00757 static int write_dcdheader(fio_fd fd, const char *remarks, int N, 
00758                     int ISTART, int NSAVC, double DELTA, int with_unitcell,
00759                     int charmm) {
00760   int out_integer;
00761   float out_float;
00762   char title_string[200];
00763   time_t cur_time;
00764   struct tm *tmbuf;
00765   char time_str[81];
00766 
00767   out_integer = 84;
00768   WRITE(fd, (char *) & out_integer, sizeof(int));
00769   strcpy(title_string, "CORD");
00770   WRITE(fd, title_string, 4);
00771   fio_write_int32(fd, 0);      /* Number of frames in file, none written yet   */
00772   fio_write_int32(fd, ISTART); /* Starting timestep                            */
00773   fio_write_int32(fd, NSAVC);  /* Timesteps between frames written to the file */
00774   fio_write_int32(fd, 0);      /* Number of timesteps in simulation            */
00775   fio_write_int32(fd, 0);      /* NAMD writes NSTEP or ISTART - NSAVC here?    */
00776   fio_write_int32(fd, 0);
00777   fio_write_int32(fd, 0);
00778   fio_write_int32(fd, 0);
00779   fio_write_int32(fd, 0);
00780   if (charmm) {
00781     out_float = DELTA;
00782     WRITE(fd, (char *) &out_float, sizeof(float));
00783     if (with_unitcell) {
00784       fio_write_int32(fd, 1);
00785     } else {
00786       fio_write_int32(fd, 0);
00787     }
00788   } else {
00789     WRITE(fd, (char *) &DELTA, sizeof(double));
00790   }
00791   fio_write_int32(fd, 0);
00792   fio_write_int32(fd, 0);
00793   fio_write_int32(fd, 0);
00794   fio_write_int32(fd, 0);
00795   fio_write_int32(fd, 0);
00796   fio_write_int32(fd, 0);
00797   fio_write_int32(fd, 0);
00798   fio_write_int32(fd, 0);
00799   if (charmm) {
00800     fio_write_int32(fd, 24); /* Pretend to be CHARMM version 24 */
00801   } else {
00802     fio_write_int32(fd, 0);
00803   }
00804   fio_write_int32(fd, 84);
00805   fio_write_int32(fd, 164);
00806   fio_write_int32(fd, 2);
00807 
00808   strncpy(title_string, remarks, 80);
00809   title_string[79] = '\0';
00810   WRITE(fd, title_string, 80);
00811 
00812   cur_time=time(NULL);
00813   tmbuf=localtime(&cur_time);
00814   strftime(time_str, 80, "REMARKS Created %d %B, %Y at %R", tmbuf);
00815   WRITE(fd, time_str, 80);
00816 
00817   fio_write_int32(fd, 164);
00818   fio_write_int32(fd, 4);
00819   fio_write_int32(fd, N);
00820   fio_write_int32(fd, 4);
00821 
00822   return DCD_SUCCESS;
00823 }
00824 
00825 
00826 /*
00827  * clean up dcd data
00828  * Input: nfixed, freeind - elements as returned by read_dcdheader
00829  * Output: None
00830  * Side effects: Space pointed to by freeind is freed if necessary.
00831  */
00832 static void close_dcd_read(int *indexes, float *fixedcoords) {
00833   free(indexes);
00834   free(fixedcoords);
00835 }
00836 
00837 
00838 
00839 
00840 
00841 static void *open_dcd_read(const char *path, const char *filetype, 
00842     int *natoms) {
00843   dcdhandle *dcd;
00844   fio_fd fd;
00845   int rc;
00846   struct stat stbuf;
00847 
00848   if (!path) return NULL;
00849 
00850   /* See if the file exists, and get its size */
00851   memset(&stbuf, 0, sizeof(struct stat));
00852   if (stat(path, &stbuf)) {
00853     printf("dcdplugin) Could not access file '%s'.\n", path);
00854     return NULL;
00855   }
00856 
00857   if (fio_open(path, FIO_READ, &fd) < 0) {
00858     printf("dcdplugin) Could not open file '%s' for reading.\n", path);
00859     return NULL;
00860   }
00861 
00862   dcd = (dcdhandle *)malloc(sizeof(dcdhandle));
00863   memset(dcd, 0, sizeof(dcdhandle));
00864   dcd->fd = fd;
00865 
00866   if ((rc = read_dcdheader(dcd->fd, &dcd->natoms, &dcd->nsets, &dcd->istart, 
00867          &dcd->nsavc, &dcd->delta, &dcd->nfixed, &dcd->freeind, 
00868          &dcd->fixedcoords, &dcd->reverse, &dcd->charmm))) {
00869     print_dcderror("read_dcdheader", rc);
00870     fio_fclose(dcd->fd);
00871     free(dcd);
00872     return NULL;
00873   }
00874 
00875   /*
00876    * Check that the file is big enough to really hold the number of sets
00877    * it claims to have.  Then we'll use nsets to keep track of where EOF
00878    * should be.
00879    */
00880   {
00881     fio_size_t ndims, firstframesize, framesize, extrablocksize;
00882     fio_size_t trjsize, filesize, curpos;
00883     int newnsets;
00884 
00885     extrablocksize = dcd->charmm & DCD_HAS_EXTRA_BLOCK ? 48 + 8 : 0;
00886     ndims = dcd->charmm & DCD_HAS_4DIMS ? 4 : 3;
00887     firstframesize = (dcd->natoms+2) * ndims * sizeof(float) + extrablocksize;
00888     framesize = (dcd->natoms-dcd->nfixed+2) * ndims * sizeof(float) 
00889       + extrablocksize;
00890 
00891     /* 
00892      * It's safe to use ftell, even though ftell returns a long, because the 
00893      * header size is < 4GB.
00894      */
00895 
00896     curpos = fio_ftell(dcd->fd); /* save current offset (end of header) */
00897 
00898 #if defined(_MSC_VER) && defined(FASTIO_NATIVEWIN32)
00899     /* the stat() call is not 64-bit savvy on Windows             */
00900     /* so we have to use the fastio fseek/ftell routines for this */
00901     /* until we add a portable filesize routine for this purpose  */
00902     fio_fseek(dcd->fd, 0, FIO_SEEK_END);       /* seek to end of file */
00903     filesize = fio_ftell(dcd->fd);
00904     fio_fseek(dcd->fd, curpos, FIO_SEEK_SET);  /* return to end of header */
00905 #else
00906     filesize = stbuf.st_size; /* this works ok on Unix machines */
00907 #endif
00908     trjsize = filesize - curpos - firstframesize;
00909     if (trjsize < 0) {
00910       printf("dcdplugin) file '%s' appears to contain no timesteps.\n", path);
00911       fio_fclose(dcd->fd);
00912       free(dcd);
00913       return NULL;
00914     }
00915 
00916     newnsets = trjsize / framesize + 1;
00917 
00918     if (dcd->nsets > 0 && newnsets != dcd->nsets) {
00919       printf("dcdplugin) Warning: DCD header claims %d frames, file size indicates there are actually %d frames\n", dcd->nsets, newnsets);
00920     }
00921 
00922     dcd->nsets = newnsets; 
00923     dcd->setsread = 0;
00924   }
00925 
00926   dcd->first = 1;
00927   dcd->x = (float *)malloc(dcd->natoms * sizeof(float));
00928   dcd->y = (float *)malloc(dcd->natoms * sizeof(float));
00929   dcd->z = (float *)malloc(dcd->natoms * sizeof(float));
00930   if (!dcd->x || !dcd->y || !dcd->z) {
00931     printf("dcdplugin) Unable to allocate space for %d atoms.\n", dcd->natoms);
00932     if (dcd->x)
00933       free(dcd->x);
00934     if (dcd->y)
00935       free(dcd->y);
00936     if (dcd->z)
00937       free(dcd->z);
00938     fio_fclose(dcd->fd);
00939     free(dcd);
00940     return NULL;
00941   }
00942   *natoms = dcd->natoms;
00943   return dcd;
00944 }
00945 
00946 
00947 static int read_next_timestep(void *v, int natoms, molfile_timestep_t *ts) {
00948   dcdhandle *dcd;
00949   int i, j, rc;
00950   float unitcell[6];
00951   unitcell[0] = unitcell[2] = unitcell[5] = 1.0f;
00952   unitcell[1] = unitcell[3] = unitcell[4] = 90.0f;
00953   dcd = (dcdhandle *)v;
00954 
00955   /* Check for EOF here; that way all EOF's encountered later must be errors */
00956   if (dcd->setsread == dcd->nsets) return MOLFILE_EOF;
00957   dcd->setsread++;
00958   if (!ts) {
00959     if (dcd->first && dcd->nfixed) {
00960       /* We can't just skip it because we need the fixed atom coordinates */
00961       rc = read_dcdstep(dcd->fd, dcd->natoms, dcd->x, dcd->y, dcd->z, 
00962           unitcell, dcd->nfixed, dcd->first, dcd->freeind, dcd->fixedcoords, 
00963              dcd->reverse, dcd->charmm);
00964       dcd->first = 0;
00965       return rc; /* XXX this needs to be updated */
00966     }
00967     dcd->first = 0;
00968     /* XXX this needs to be changed */
00969     return skip_dcdstep(dcd->fd, dcd->natoms, dcd->nfixed, dcd->charmm);
00970   }
00971   rc = read_dcdstep(dcd->fd, dcd->natoms, dcd->x, dcd->y, dcd->z, unitcell,
00972              dcd->nfixed, dcd->first, dcd->freeind, dcd->fixedcoords, 
00973              dcd->reverse, dcd->charmm);
00974   dcd->first = 0;
00975   if (rc < 0) {  
00976     print_dcderror("read_dcdstep", rc);
00977     return MOLFILE_ERROR;
00978   }
00979 
00980   /* copy timestep data from plugin-local buffers to VMD's buffer */
00981   /* XXX 
00982    *   This code is still the root of all evil.  Just doing this extra copy
00983    *   cuts the I/O rate of the DCD reader from 728 MB/sec down to
00984    *   394 MB/sec when reading from a ram filesystem.  
00985    *   For a physical disk filesystem, the I/O rate goes from 
00986    *   187 MB/sec down to 122 MB/sec.  Clearly this extra copy has to go.
00987    */
00988   {
00989     int natoms = dcd->natoms;
00990     float *nts = ts->coords;
00991     const float *bufx = dcd->x;
00992     const float *bufy = dcd->y;
00993     const float *bufz = dcd->z;
00994 
00995     for (i=0, j=0; i<natoms; i++, j+=3) {
00996       nts[j    ] = bufx[i];
00997       nts[j + 1] = bufy[i];
00998       nts[j + 2] = bufz[i];
00999     }
01000   }
01001 
01002   ts->A = unitcell[0];
01003   ts->B = unitcell[2];
01004   ts->C = unitcell[5];
01005 
01006   if (unitcell[1] >= -1.0 && unitcell[1] <= 1.0 &&
01007       unitcell[3] >= -1.0 && unitcell[3] <= 1.0 &&
01008       unitcell[4] >= -1.0 && unitcell[4] <= 1.0) {
01009     /* This file was generated by CHARMM, or by NAMD > 2.5, with the angle */
01010     /* cosines of the periodic cell angles written to the DCD file.        */ 
01011     /* This formulation improves rounding behavior for orthogonal cells    */
01012     /* so that the angles end up at precisely 90 degrees, unlike acos().   */
01013     ts->alpha = 90.0 - asin(unitcell[4]) * 90.0 / M_PI_2; /* cosBC */
01014     ts->beta  = 90.0 - asin(unitcell[3]) * 90.0 / M_PI_2; /* cosAC */
01015     ts->gamma = 90.0 - asin(unitcell[1]) * 90.0 / M_PI_2; /* cosAB */
01016   } else {
01017     /* This file was likely generated by NAMD 2.5 and the periodic cell    */
01018     /* angles are specified in degrees rather than angle cosines.          */
01019     ts->alpha = unitcell[4]; /* angle between B and C */
01020     ts->beta  = unitcell[3]; /* angle between A and C */
01021     ts->gamma = unitcell[1]; /* angle between A and B */
01022   }
01023  
01024   return MOLFILE_SUCCESS;
01025 }
01026  
01027 
01028 static void close_file_read(void *v) {
01029   dcdhandle *dcd = (dcdhandle *)v;
01030   close_dcd_read(dcd->freeind, dcd->fixedcoords);
01031   fio_fclose(dcd->fd);
01032   free(dcd->x);
01033   free(dcd->y);
01034   free(dcd->z);
01035   free(dcd); 
01036 }
01037 
01038 
01039 static void *open_dcd_write(const char *path, const char *filetype, 
01040     int natoms) {
01041   dcdhandle *dcd;
01042   fio_fd fd;
01043   int rc;
01044   int istart, nsavc;
01045   double delta;
01046   int with_unitcell;
01047   int charmm;
01048 
01049   if (fio_open(path, FIO_WRITE, &fd) < 0) {
01050     printf("dcdplugin) Could not open file '%s' for writing\n", path);
01051     return NULL;
01052   }
01053 
01054   dcd = (dcdhandle *)malloc(sizeof(dcdhandle));
01055   memset(dcd, 0, sizeof(dcdhandle));
01056   dcd->fd = fd;
01057 
01058   istart = 0;             /* starting timestep of DCD file                  */
01059   nsavc = 1;              /* number of timesteps between written DCD frames */
01060   delta = 1.0;            /* length of a timestep                           */
01061 
01062   if (getenv("VMDDCDWRITEXPLORFORMAT") != NULL) {
01063     with_unitcell = 0;      /* no unit cell info */
01064     charmm = DCD_IS_XPLOR;  /* X-PLOR format */
01065     printf("dcdplugin) WARNING: Writing DCD file in X-PLOR format, \n");
01066     printf("dcdplugin) WARNING: unit cell information will be lost!\n");
01067   } else {
01068     with_unitcell = 1;      /* contains unit cell infor (Charmm format) */
01069     charmm = DCD_IS_CHARMM; /* charmm-formatted DCD file                */ 
01070     if (with_unitcell) 
01071       charmm |= DCD_HAS_EXTRA_BLOCK;
01072   }
01073  
01074   rc = write_dcdheader(dcd->fd, "Created by DCD plugin", natoms, 
01075                        istart, nsavc, delta, with_unitcell, charmm);
01076 
01077   if (rc < 0) {
01078     print_dcderror("write_dcdheader", rc);
01079     fio_fclose(dcd->fd);
01080     free(dcd);
01081     return NULL;
01082   }
01083 
01084   dcd->natoms = natoms;
01085   dcd->nsets = 0;
01086   dcd->istart = istart;
01087   dcd->nsavc = nsavc;
01088   dcd->with_unitcell = with_unitcell;
01089   dcd->charmm = charmm;
01090   dcd->x = (float *)malloc(natoms * sizeof(float));
01091   dcd->y = (float *)malloc(natoms * sizeof(float));
01092   dcd->z = (float *)malloc(natoms * sizeof(float));
01093   return dcd;
01094 }
01095 
01096 
01097 static int write_timestep(void *v, const molfile_timestep_t *ts) { 
01098   dcdhandle *dcd = (dcdhandle *)v;
01099   int i, rc, curstep;
01100   float *pos = ts->coords;
01101   double unitcell[6];
01102   unitcell[0] = unitcell[2] = unitcell[5] = 1.0f;
01103   unitcell[1] = unitcell[3] = unitcell[4] = 90.0f;
01104 
01105   /* copy atom coords into separate X/Y/Z arrays for writing */
01106   for (i=0; i<dcd->natoms; i++) {
01107     dcd->x[i] = *(pos++); 
01108     dcd->y[i] = *(pos++); 
01109     dcd->z[i] = *(pos++); 
01110   }
01111   dcd->nsets++;
01112   curstep = dcd->istart + dcd->nsets * dcd->nsavc;
01113 
01114   unitcell[0] = ts->A;
01115   unitcell[2] = ts->B;
01116   unitcell[5] = ts->C;
01117   unitcell[1] = sin((M_PI_2 / 90.0) * (90.0 - ts->gamma)); /* cosAB */
01118   unitcell[3] = sin((M_PI_2 / 90.0) * (90.0 - ts->beta));  /* cosAC */
01119   unitcell[4] = sin((M_PI_2 / 90.0) * (90.0 - ts->alpha)); /* cosBC */
01120 
01121   rc = write_dcdstep(dcd->fd, dcd->nsets, curstep, dcd->natoms, 
01122                      dcd->x, dcd->y, dcd->z,
01123                      dcd->with_unitcell ? unitcell : NULL,
01124                      dcd->charmm);
01125   if (rc < 0) {
01126     print_dcderror("write_dcdstep", rc);
01127     return MOLFILE_ERROR;
01128   }
01129 
01130   return MOLFILE_SUCCESS;
01131 }
01132 
01133 static void close_file_write(void *v) {
01134   dcdhandle *dcd = (dcdhandle *)v;
01135   fio_fclose(dcd->fd);
01136   free(dcd->x);
01137   free(dcd->y);
01138   free(dcd->z);
01139   free(dcd);
01140 }
01141 
01142 
01143 /*
01144  * Initialization stuff here
01145  */
01146 static molfile_plugin_t plugin;
01147 
01148 VMDPLUGIN_API int VMDPLUGIN_init() {
01149   memset(&plugin, 0, sizeof(molfile_plugin_t));
01150   plugin.abiversion = vmdplugin_ABIVERSION;
01151   plugin.type = MOLFILE_PLUGIN_TYPE;
01152   plugin.name = "dcd";
01153   plugin.prettyname = "CHARMM,NAMD,XPLOR DCD Trajectory";
01154   plugin.author = "Axel Kohlmeyer, Justin Gullingsrud, John Stone";
01155   plugin.majorv = 1;
01156   plugin.minorv = 12;
01157   plugin.is_reentrant = VMDPLUGIN_THREADSAFE;
01158   plugin.filename_extension = "dcd";
01159   plugin.open_file_read = open_dcd_read;
01160   plugin.read_next_timestep = read_next_timestep;
01161   plugin.close_file_read = close_file_read;
01162   plugin.open_file_write = open_dcd_write;
01163   plugin.write_timestep = write_timestep;
01164   plugin.close_file_write = close_file_write;
01165   return VMDPLUGIN_SUCCESS;
01166 }
01167 
01168 VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
01169   (*cb)(v, (vmdplugin_t *)&plugin);
01170   return VMDPLUGIN_SUCCESS;
01171 }
01172 
01173 VMDPLUGIN_API int VMDPLUGIN_fini() {
01174   return VMDPLUGIN_SUCCESS;
01175 }
01176 
01177   
01178 #ifdef TEST_DCDPLUGIN
01179 
01180 #include <sys/time.h>
01181 
01182 /* get the time of day from the system clock, and store it (in seconds) */
01183 double time_of_day(void) {
01184 #if defined(_MSC_VER)
01185   double t;
01186 
01187   t = GetTickCount();
01188   t = t / 1000.0;
01189 
01190   return t;
01191 #else
01192   struct timeval tm;
01193   struct timezone tz;
01194 
01195   gettimeofday(&tm, &tz);
01196   return((double)(tm.tv_sec) + (double)(tm.tv_usec)/1000000.0);
01197 #endif
01198 }
01199 
01200 int main(int argc, char *argv[]) {
01201   molfile_timestep_t timestep;
01202   void *v;
01203   dcdhandle *dcd;
01204   int i, natoms;
01205   float sizeMB =0.0, totalMB = 0.0;
01206   double starttime, endtime, totaltime = 0.0;
01207 
01208   while (--argc) {
01209     ++argv; 
01210     natoms = 0;
01211     v = open_dcd_read(*argv, "dcd", &natoms);
01212     if (!v) {
01213       fprintf(stderr, "main) open_dcd_read failed for file %s\n", *argv);
01214       return 1;
01215     }
01216     dcd = (dcdhandle *)v;
01217     sizeMB = ((natoms * 3.0) * dcd->nsets * 4.0) / (1024.0 * 1024.0);
01218     totalMB += sizeMB; 
01219     printf("main) file: %s\n", *argv);
01220     printf("  %d atoms, %d frames, size: %6.1fMB\n", natoms, dcd->nsets, sizeMB);
01221 
01222     starttime = time_of_day();
01223     timestep.coords = (float *)malloc(3*sizeof(float)*natoms);
01224     for (i=0; i<dcd->nsets; i++) {
01225       int rc = read_next_timestep(v, natoms, &timestep);
01226       if (rc) {
01227         fprintf(stderr, "error in read_next_timestep on frame %d\n", i);
01228         return 1;
01229       }
01230     }
01231     endtime = time_of_day();
01232     close_file_read(v);
01233     totaltime += endtime - starttime;
01234     printf("  Time: %5.1f seconds\n", endtime - starttime);
01235     printf("  Speed: %5.1f MB/sec, %5.1f timesteps/sec\n", sizeMB / (endtime - starttime), (dcd->nsets / (endtime - starttime)));
01236   }
01237   printf("Overall Size: %6.1f MB\n", totalMB);
01238   printf("Overall Time: %6.1f seconds\n", totaltime);
01239   printf("Overall Speed: %5.1f MB/sec\n", totalMB / totaltime);
01240   return 0;
01241 }
01242       
01243 #endif
01244 

Generated on Thu Apr 24 02:34:47 2014 for VMD Plugins (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002