dcdlib.C

Go to the documentation of this file.
00001 
00007 /*
00008    dcdlib contains C routines for reading and writing binary DCD
00009    files.  The output format of these files is based on binary FORTRAN
00010    output, so its pretty ugly.  If you are squimish, don't look!
00011 */
00012 
00013 #include "dcdlib.h"
00014 
00015 #ifndef OUTPUT_SINGLE_FILE
00016 #define OUTPUT_SINGLE_FILE 1
00017 // including Output.h causes compile error on BlueGeneQ
00018 #endif
00019 
00020 #ifdef WIN32
00021 #define access(PATH,MODE) _access(PATH,00)
00022 #endif
00023 
00024 #define NAMD_write NAMD_write64
00025 // same as write, only does error checking internally
00026 void NAMD_write(int fd, const char *buf, size_t count) {
00027   while ( count ) {
00028 #if defined(WIN32) && !defined(__CYGWIN__)
00029     long retval = _write(fd,buf,count);
00030 #else
00031     ssize_t retval = write(fd,buf,count);
00032 #endif
00033     if ( retval < 0 && errno == EINTR ) retval = 0;
00034     if ( retval < 0 ) NAMD_die(strerror(errno));
00035     if ( retval > count ) NAMD_bug("extra bytes written in NAMD_write64()");
00036     buf += retval;
00037     count -= retval;
00038   }
00039 }
00040 
00041 #ifdef WIN32
00042 #define LSEEK _lseeki64
00043 #define READ _read
00044 #else
00045 #define LSEEK lseek
00046 #define READ read
00047 #endif
00048 
00049 OFF_T NAMD_seek(int file, OFF_T offset, int whence) {
00050   OFF_T retval = LSEEK(file, offset, whence);
00051   if ( retval < 0 ) NAMD_err("seek failed while writing DCD file");
00052   if ( whence == SEEK_SET && retval != offset ) {
00053     char buf[256];
00054     sprintf(buf, "seek failed while writing DCD file: SEEK_SET %lld returned %lld\n", offset, retval);
00055     NAMD_die(buf);
00056   }
00057   return retval;
00058 }
00059 
00060 #undef LSEEK
00061 #define LSEEK NAMD_seek
00062 
00063 #ifndef O_LARGEFILE
00064 #define O_LARGEFILE 0x0
00065 #endif
00066 
00067 
00068 /************************************************************************/
00069 /*                                                                      */
00070 /*                              MACRO CHECK_FREAD                       */
00071 /*                                                                      */
00072 /*      CHECK_FREAD is used to check the status of a file read.  It     */
00073 /*   is passed the return code from read and a string to print out if   */
00074 /*   an error is detected.  If an error is found, an error message is   */
00075 /*   printed out and the program terminates.  This was made into a      */
00076 /*   macro because it had to be done over and over and over . . .       */
00077 /*                                                                      */
00078 /************************************************************************/
00079 
00080 #define CHECK_FREAD(X, msg)  if (X==-1) \
00081                              { \
00082                                 return(DCD_BADREAD); \
00083                              }
00084 
00085 /************************************************************************/
00086 /*                                                                      */
00087 /*                      MACRO CHECK_FEOF                                */
00088 /*                                                                      */
00089 /*      CHECK_FEOF checks for an abnormal EOF on a file read.  It is    */
00090 /*   passed the return status from read and a message to print on error.*/
00091 /*   if an EOF is detected, the error message is printed and the program*/
00092 /*   terminates.  This is used for reads that should find something in  */
00093 /*   the file.                                                          */
00094 /*                                                                      */
00095 /************************************************************************/
00096 
00097 #define CHECK_FEOF(X, msg)  if (X==0) \
00098                              { \
00099                                 return(DCD_BADEOF); \
00100                              }
00101 
00102 
00103 /*                      GLOBAL VARIABLES                                */
00104 
00105 void pad(char *s, int len)
00106 {
00107         int curlen;
00108         int i;
00109 
00110         curlen=strlen(s);
00111 
00112         if (curlen>len)
00113         {
00114                 s[len]='\0';
00115                 return;
00116         }
00117 
00118         for (i=curlen; i<len; i++)
00119         {
00120                 s[i]=' ';
00121         }
00122 
00123         s[i]='\0';
00124 }
00125 
00126 /************************************************************************/
00127 /*                                                                      */
00128 /*                      FUNCTION open_dcd_read                          */
00129 /*                                                                      */
00130 /*   INPUTS:                                                            */
00131 /*      filename - filename to read in                                  */
00132 /*                                                                      */
00133 /*   OUTPUTS:                                                           */
00134 /*      an open filedesriptor (integer) is returned if the call is      */
00135 /*   successful, otherwise a negative number is returned                */
00136 /*                                                                      */
00137 /*      open_dcd_read opens a dcd file for input.  It first does a check*/
00138 /*   to see if the file really exists.  If the file does not exist,     */
00139 /*   a DCD_DNE is returned, if the file exists but can' be opened for   */
00140 /*   some reason, a DCD_OPENFAILED is returned.  If the open is         */
00141 /*   successful, the file descriptor is returned.                       */
00142 /*                                                                      */
00143 /************************************************************************/
00144 
00145 int open_dcd_read(char *filename)
00146 
00147 {
00148         int dcdfd;              /*  file descriptor for dcd file        */
00149 
00150         /*  Do a stat just to see if the file really exists     */
00151         if (access(filename, F_OK) != 0)
00152         {
00153                 if (errno == ENOENT)
00154                 {
00155                         return(DCD_DNE);
00156                 }
00157         }
00158 
00159         /*  Try and open the file                               */
00160 #ifdef WIN32
00161         dcdfd=_open(filename, O_RDONLY|O_BINARY|O_LARGEFILE);
00162 #else
00163         dcdfd=open(filename, O_RDONLY|O_LARGEFILE);
00164 #endif
00165 
00166         if (dcdfd == -1)
00167         {
00168                 return(DCD_OPENFAILED);
00169         }
00170 
00171         return(dcdfd);
00172 }
00173 
00174 /****************************************************************/
00175 /*                                                              */
00176 /*                      FUNCTION read_dcdheader                 */
00177 /*                                                              */
00178 /*   INPUTS:                                                    */
00179 /*      fd - file descriptor for the dcd file                   */
00180 /*      N - Number of atoms                                     */
00181 /*      NSET - Number of sets of coordinates                    */
00182 /*      ISTART - Starting timestep of DCD file                  */
00183 /*      NSAVC - Timesteps between DCD saves                     */
00184 /*      DELTA - length of a timestep                            */
00185 /*                                                              */
00186 /*   OUTPUTS:                                                   */
00187 /*      N, NSET, ISTART, NSAVC, and DELTA are returned populated*/
00188 /*   a 0 is returned if the call is successful, or a negative   */
00189 /*   number if errors are detected                              */
00190 /*                                                              */
00191 /*      read_header reads in the header from a DCD file and     */
00192 /*   returns the timestep size and the number of atoms in the   */
00193 /*   system.  A 0 is returned if the header is successfully     */
00194 /*   read, or a negative error code is returned otherwise.      */
00195 /*                                                              */
00196 /****************************************************************/
00197 
00198 #if 0
00199 int read_dcdheader(int fd, int *N, int *NSET, int *ISTART, 
00200                    int *NSAVC, double *DELTA, int *NAMNF, 
00201                    int **FREEINDEXES)
00202 
00203 {
00204         int32 input_integer;    /*  Integer buffer space        */
00205         float input_float;
00206         char bigbuf[256];       /*  A large string buffer       */
00207         int ret_val;            /*  Return value from read      */
00208         int i;                  /*  Loop counter                */
00209         char HDR[5];            /*  Header = "CORD"             */
00210         int32 I;
00211         int32 NTITLE;
00212 
00213         /*  First thing in the file should be an 84             */
00214         ret_val = read(fd, &input_integer, sizeof(int32));
00215 
00216         CHECK_FREAD(ret_val, "reading first int from dcd file");
00217         CHECK_FEOF(ret_val, "reading first int from dcd file");
00218 
00219         if (input_integer != 84)
00220         {
00221                 return(DCD_BADFORMAT);
00222         }
00223 
00224         /*  Read in the string "CORD"                   */
00225         ret_val = read(fd, HDR, 4);
00226 
00227         CHECK_FREAD(ret_val, "reading CORD from dcd file");
00228         CHECK_FEOF(ret_val, "reading CORD from dcd file");
00229 
00230         HDR[ret_val]='\0';
00231 
00232         if (strcmp(HDR, "CORD") != 0)
00233         {
00234                 return(DCD_BADFORMAT);
00235         }
00236 
00237         /*  Read in the number of Sets of coordinates, NSET  */
00238         ret_val = read(fd, &input_integer, sizeof(int32));
00239         *NSET = input_integer;
00240 
00241         CHECK_FREAD(ret_val, "reading NSET");
00242         CHECK_FEOF(ret_val, "reading NSET");
00243 
00244         /*  Read in ISTART, the starting timestep            */
00245         ret_val = read(fd, &input_integer, sizeof(int32));
00246         *ISTART = input_integer;
00247 
00248         CHECK_FREAD(ret_val, "reading ISTART");
00249         CHECK_FEOF(ret_val, "reading ISTART");
00250 
00251         /*  Read in NSAVC, the number of timesteps between   */
00252         /*  dcd saves                                        */
00253         ret_val = read(fd, &input_integer, sizeof(int32));
00254         *NSAVC = input_integer;
00255 
00256         CHECK_FREAD(ret_val, "reading NSAVC");
00257         CHECK_FEOF(ret_val, "reading NSAVC");
00258 
00259 
00260         /*  Skip blank integers                              */
00261         for (i=0; i<5; i++)
00262         {
00263                 ret_val = read(fd, &input_integer, sizeof(int32));
00264 
00265                 CHECK_FREAD(ret_val, "reading I");
00266                 CHECK_FEOF(ret_val, "reading I");
00267         }
00268 
00269         /*  Read NAMNF, the number of free atoms             */
00270         ret_val = read(fd, &input_integer, sizeof(int32));
00271         *NAMNF = input_integer;
00272 
00273         CHECK_FREAD(ret_val, "reading NAMNF");
00274         CHECK_FEOF(ret_val, "reading NAMNF");
00275 
00276         /*  Read in the timestep, DELTA                         */
00277         ret_val = read(fd, &input_float, sizeof(float));
00278 
00279         CHECK_FREAD(ret_val, "reading DELTA");
00280         CHECK_FEOF(ret_val, "reading DELTA");
00281         *DELTA = input_float;
00282 
00283         /*  Skip blank integers                                 */
00284         for (i=0; i<10; i++)
00285         {
00286                 ret_val = read(fd, &I, sizeof(int32));
00287 
00288                 CHECK_FREAD(ret_val, "reading I");
00289                 CHECK_FEOF(ret_val, "reading I");
00290         }
00291 
00292         /*  Get the end size of the first block                 */
00293         ret_val = read(fd, &input_integer, sizeof(int32));
00294 
00295         CHECK_FREAD(ret_val, "reading second 84 from dcd file");
00296         CHECK_FEOF(ret_val, "reading second 84 from dcd file");
00297 
00298         if (input_integer != 84)
00299         {
00300                 return(DCD_BADFORMAT);
00301         }
00302 
00303         /*  Read in the size of the next block                  */
00304         ret_val = read(fd, &input_integer, sizeof(int32));
00305 
00306         CHECK_FREAD(ret_val, "reading size of title block");
00307         CHECK_FEOF(ret_val, "reading size of title block");
00308 
00309         if ( ((input_integer-4)%80) == 0)
00310         {
00311                 /*  Read NTITLE, the number of 80 characeter    */
00312                 /*  title strings there are                     */
00313                 ret_val = read(fd, &NTITLE, sizeof(int32));
00314 
00315                 CHECK_FREAD(ret_val, "reading NTITLE");
00316                 CHECK_FEOF(ret_val, "reading NTITLE");
00317 
00318                 for (i=0; i<NTITLE; i++)
00319                 {
00320                         ret_val = read(fd, bigbuf, 80);
00321 
00322                         CHECK_FREAD(ret_val, "reading TITLE");
00323                         CHECK_FEOF(ret_val, "reading TITLE");
00324                 }
00325 
00326                 /*  Get the ending size for this block          */
00327                 ret_val = read(fd, &input_integer, sizeof(int32));
00328 
00329                 CHECK_FREAD(ret_val, "reading size of title block");
00330                 CHECK_FEOF(ret_val, "reading size of title block");
00331         }
00332         else
00333         {
00334                 return(DCD_BADFORMAT);
00335         }
00336 
00337         /*  Read in an 4                                */
00338         ret_val = read(fd, &input_integer, sizeof(int32));
00339 
00340         CHECK_FREAD(ret_val, "reading an 4");
00341         CHECK_FEOF(ret_val, "reading an 4");
00342 
00343         if (input_integer != 4)
00344         {
00345                 return(DCD_BADFORMAT);
00346         }
00347 
00348         /*  Read in the number of atoms                 */
00349         ret_val = read(fd, &input_integer, sizeof(int32));
00350         *N = input_integer;
00351 
00352         CHECK_FREAD(ret_val, "reading number of atoms");
00353         CHECK_FEOF(ret_val, "reading number of atoms");
00354 
00355         /*  Read in an 4                                */
00356         ret_val = read(fd, &input_integer, sizeof(int32));
00357 
00358         CHECK_FREAD(ret_val, "reading an 4");
00359         CHECK_FEOF(ret_val, "reading an 4");
00360 
00361         if (input_integer != 4)
00362         {
00363                 return(DCD_BADFORMAT);
00364         }
00365 
00366         if (*NAMNF != 0)
00367         {
00368                 (*FREEINDEXES) = new int[(*N)-(*NAMNF)];
00369                 int32 *freeindexes32 = new int32[(*N)-(*NAMNF)];
00370 
00371                 if (*FREEINDEXES == NULL || freeindexes32 == NULL)
00372                         return(DCD_BADMALLOC);
00373         
00374                 /*  Read in an size                             */
00375                 ret_val = read(fd, &input_integer, sizeof(int32));
00376 
00377                 CHECK_FREAD(ret_val, "reading size of index array");
00378                 CHECK_FEOF(ret_val, "reading size of index array");
00379 
00380                 if (input_integer != ((*N)-(*NAMNF))*4)
00381                 {
00382                         return(DCD_BADFORMAT);
00383                 }
00384                 
00385                 ret_val = read(fd, freeindexes32, ((*N)-(*NAMNF))*sizeof(int32));
00386 
00387                 CHECK_FREAD(ret_val, "reading size of index array");
00388                 CHECK_FEOF(ret_val, "reading size of index array");
00389 
00390                 for (i=0; i<((*N)-(*NAMNF)); ++i)
00391                         (*FREEINDEXES)[i] = freeindexes32[i];
00392 
00393                 delete [] freeindexes32;
00394 
00395                 ret_val = read(fd, &input_integer, sizeof(int32));
00396 
00397                 CHECK_FREAD(ret_val, "reading size of index array");
00398                 CHECK_FEOF(ret_val, "reading size of index array");
00399 
00400                 if (input_integer != ((*N)-(*NAMNF))*4)
00401                 {
00402                         return(DCD_BADFORMAT);
00403                 }
00404         }
00405 
00406         return(0);
00407 }
00408 
00409 /************************************************************************/
00410 /*                                                                      */
00411 /*                      FUNCTION read_dcdstep                           */
00412 /*                                                                      */
00413 /*   INPUTS:                                                            */
00414 /*      fd - file descriptor to use                                     */
00415 /*      N - Number of atoms                                             */
00416 /*      X - array of X coordinates                                      */
00417 /*      Y - array of Y coordinates                                      */
00418 /*      Z - array of Z coordinates                                      */
00419 /*      num_fixed - Number of fixed atoms                               */
00420 /*      first - flag 1->first time called                               */
00421 /*      indexes - indexes for free atoms                                */
00422 /*                                                                      */
00423 /*   OUTPUTS:                                                           */
00424 /*      read step populates the arrays X, Y, Z and returns a 0 if the   */
00425 /*   read is successful.  If the read fails, a negative error code is   */
00426 /*   returned.                                                          */
00427 /*                                                                      */
00428 /*      read_step reads in the coordinates from one step.  It places    */
00429 /*   these coordinates into the arrays X, Y, and Z.                     */
00430 /*                                                                      */
00431 /************************************************************************/
00432 
00433 int read_dcdstep(int fd, int N, float *X, float *Y, float *Z, int num_fixed,
00434                  int first, int *indexes)
00435 
00436 {
00437         int ret_val;            /*  Return value from read              */
00438         int32 input_integer;    /*  Integer buffer space                */
00439         int i;                  /*  Loop counter                        */
00440         static float *tmpX;
00441 
00442         if (first && num_fixed)
00443         {
00444                 tmpX = new float[N-num_fixed];
00445 
00446                 if (tmpX==NULL)
00447                 {
00448                         return(DCD_BADMALLOC);
00449                 }
00450         }
00451 
00452         /*  Get the first size from the file                            */
00453         ret_val = read(fd, &input_integer, sizeof(int32));
00454 
00455         CHECK_FREAD(ret_val, "reading number of atoms at begining of step");
00456 
00457         /*  See if we've reached the end of the file                    */
00458         if (ret_val == 0)
00459         {
00460                 delete [] tmpX;
00461 
00462                 return(-1);
00463         }
00464 
00465         if ( (num_fixed==0) || first)
00466         {
00467                 if (input_integer != 4*N)
00468                 {
00469                         return(DCD_BADFORMAT);
00470                 }
00471 
00472                 ret_val = read(fd, X, 4*N);
00473         
00474                 CHECK_FREAD(ret_val, "reading X array");
00475                 CHECK_FEOF(ret_val, "reading X array");
00476 
00477                 ret_val = read(fd, &input_integer, sizeof(int32));
00478 
00479                 CHECK_FREAD(ret_val, "reading number of atoms after X array");
00480                 CHECK_FEOF(ret_val, "reading number of atoms after X array");
00481 
00482                 if (input_integer != 4*N)
00483                 {
00484                         return(DCD_BADFORMAT);
00485                 }
00486 
00487                 ret_val = read(fd, &input_integer, sizeof(int32));
00488 
00489                 CHECK_FREAD(ret_val, "reading number of atoms after X array");
00490                 CHECK_FEOF(ret_val, "reading number of atoms after X array");
00491 
00492                 if (input_integer != 4*N)
00493                 {
00494                         return(DCD_BADFORMAT);
00495                 }
00496 
00497                 ret_val = read(fd, Y, 4*N);
00498 
00499                 CHECK_FREAD(ret_val, "reading Y array");
00500                 CHECK_FEOF(ret_val, "reading Y array");
00501 
00502                 ret_val = read(fd, &input_integer, sizeof(int32));
00503 
00504                 CHECK_FREAD(ret_val, "reading number of atoms after Y array");
00505                 CHECK_FEOF(ret_val, "reading number of atoms after Y array");
00506 
00507                 if (input_integer != 4*N)
00508                 {
00509                         return(DCD_BADFORMAT);
00510                 }
00511 
00512                 ret_val = read(fd, &input_integer, sizeof(int32));
00513 
00514                 CHECK_FREAD(ret_val, "reading number of atoms after Y array");
00515                 CHECK_FEOF(ret_val, "reading number of atoms after Y array");
00516 
00517                 if (input_integer != 4*N)
00518                 {
00519                         return(DCD_BADFORMAT);
00520                 }
00521 
00522                 ret_val = read(fd, Z, 4*N);
00523 
00524                 CHECK_FREAD(ret_val, "reading Z array");
00525                 CHECK_FEOF(ret_val, "reading Z array");
00526 
00527                 ret_val = read(fd, &input_integer, sizeof(int32));
00528 
00529                 CHECK_FREAD(ret_val, "reading number of atoms after Z array");
00530                 CHECK_FEOF(ret_val, "reading number of atoms after Z array");
00531 
00532                 if (input_integer != 4*N)
00533                 {
00534                         return(DCD_BADFORMAT);
00535                 }
00536         }
00537         else
00538         {
00539                 if (input_integer != 4*(N-num_fixed))
00540                 {
00541                         return(DCD_BADFORMAT);
00542                 }
00543 
00544                 ret_val = read(fd, tmpX, 4*(N-num_fixed));
00545         
00546                 CHECK_FREAD(ret_val, "reading Xtmp array");
00547                 CHECK_FEOF(ret_val, "reading Xtmp array");
00548 
00549                 for (i=0; i<N-num_fixed; i++)
00550                 {
00551                         X[indexes[i]-1]=tmpX[i];
00552                 }
00553 
00554                 ret_val = read(fd, &input_integer, sizeof(int32));
00555 
00556                 CHECK_FREAD(ret_val, "reading number of atoms after X array");
00557                 CHECK_FEOF(ret_val, "reading number of atoms after X array");
00558 
00559                 if (input_integer != 4*(N-num_fixed))
00560                 {
00561                         return(DCD_BADFORMAT);
00562                 }
00563 
00564                 ret_val = read(fd, &input_integer, sizeof(int32));
00565 
00566                 if (input_integer != 4*(N-num_fixed))
00567                 {
00568                         return(DCD_BADFORMAT);
00569                 }
00570 
00571                 ret_val = read(fd, tmpX, 4*(N-num_fixed));
00572         
00573                 CHECK_FREAD(ret_val, "reading Xtmp array");
00574                 CHECK_FEOF(ret_val, "reading Xtmp array");
00575 
00576                 for (i=0; i<N-num_fixed; i++)
00577                 {
00578                         Y[indexes[i]-1]=tmpX[i];
00579                 }
00580 
00581                 ret_val = read(fd, &input_integer, sizeof(int32));
00582 
00583                 CHECK_FREAD(ret_val, "reading number of atoms after Y array");
00584                 CHECK_FEOF(ret_val, "reading number of atoms after Y array");
00585 
00586                 if (input_integer != 4*(N-num_fixed))
00587                 {
00588                         return(DCD_BADFORMAT);
00589                 }
00590 
00591                 ret_val = read(fd, &input_integer, sizeof(int32));
00592 
00593                 if (input_integer != 4*(N-num_fixed))
00594                 {
00595                         return(DCD_BADFORMAT);
00596                 }
00597 
00598                 ret_val = read(fd, tmpX, 4*(N-num_fixed));
00599         
00600                 CHECK_FREAD(ret_val, "reading Xtmp array");
00601                 CHECK_FEOF(ret_val, "reading Xtmp array");
00602 
00603                 for (i=0; i<N-num_fixed; i++)
00604                 {
00605                         Z[indexes[i]-1]=tmpX[i];
00606                 }
00607 
00608                 ret_val = read(fd, &input_integer, sizeof(int32));
00609 
00610                 CHECK_FREAD(ret_val, "reading number of atoms after Z array");
00611                 CHECK_FEOF(ret_val, "reading number of atoms after Z array");
00612 
00613                 if (input_integer != 4*(N-num_fixed))
00614                 {
00615                         return(DCD_BADFORMAT);
00616                 }
00617         }
00618 
00619         return(0);
00620 }
00621 #endif
00622 
00623 #ifndef OUTPUT_SINGLE_FILE
00624 #error OUTPUT_SINGLE_FILE not defined!
00625 #endif
00626 
00627 #if OUTPUT_SINGLE_FILE
00628 #define NFILE_POS  ((OFF_T) 8)
00629 #define NPRIV_POS  ((OFF_T) 12)
00630 #define NSAVC_POS  ((OFF_T) 16)
00631 #define NSTEP_POS  ((OFF_T) 20)
00632 #else
00633 //Need to consider extra fields: 
00634 //1. magic number (int32)
00635 //2. file version (float)
00636 //3. number of files that contain portion of data (int32)
00637 //So the total offset is 12 bytes.
00638 #define NFILE_POS  ((OFF_T) 20)
00639 #define NPRIV_POS  ((OFF_T) 24)
00640 #define NSAVC_POS  ((OFF_T) 28)
00641 #define NSTEP_POS  ((OFF_T) 32)
00642 #endif
00643 
00644 /*********************************************************************/
00645 /*                                                                   */
00646 /*                      FUNCTION open_dcd_write                      */
00647 /*                                                                   */
00648 /*   INPUTS:                                                         */
00649 /*      dcdfile - Name of the dcd file                               */
00650 /*                                                                   */
00651 /*   OUTPUTS:                                                        */
00652 /*      returns an open file descriptor for writing                  */
00653 /*                                                                   */
00654 /*      This function will open a dcd file for writing.  It takes    */
00655 /*   the filename to open as its only argument.  It will return a    */
00656 /*   valid file descriptor if successful or DCD_OPENFAILED if the    */
00657 /*   open fails for some reason.  If the file specifed already       */
00658 /*   exists, it is renamed by appending .BAK to it.                  */
00659 /*                                                                   */
00660 /*********************************************************************/
00661 
00662 int open_dcd_write(const char *dcdname)
00663 
00664 {
00665         int dcdfd;
00666         NAMD_backup_file(dcdname,".BAK");
00667 
00668 #ifdef WIN32
00669         while ( (dcdfd = _open(dcdname, O_RDWR|O_CREAT|O_EXCL|O_BINARY|O_LARGEFILE,
00670                                 _S_IREAD|_S_IWRITE)) < 0)
00671 #else
00672 #ifdef NAMD_NO_O_EXCL
00673         while ( (dcdfd = open(dcdname, O_RDWR|O_CREAT|O_TRUNC|O_LARGEFILE,
00674 #else
00675         while ( (dcdfd = open(dcdname, O_RDWR|O_CREAT|O_EXCL|O_LARGEFILE,
00676 #endif
00677                                 S_IRUSR|S_IWUSR|S_IRGRP|S_IROTH)) < 0)
00678 #endif
00679         {
00680                 if ( errno != EINTR ) return(DCD_OPENFAILED);
00681         }
00682 
00683         return(dcdfd);
00684 }
00685 
00686 // non master nodes open an existing file
00687 // master will have taken care of backup and creation before the
00688 // broadcast which triggers the slaves
00689 int open_dcd_write_par_slave(char *dcdname)
00690 
00691 {
00692 #if OUTPUT_SINGLE_FILE
00693         //In the case of single file, the dcd output by slaves has been created
00694         //by the master, so the dcd file doesn't need to be created again and
00695         //backed up. --Chao Mei
00696         int dcdfd;
00697 #ifdef WIN32
00698         while ( (dcdfd = _open(dcdname, O_WRONLY|O_BINARY|O_LARGEFILE,
00699                                 _S_IREAD|_S_IWRITE)) < 0)
00700 #else
00701         while ( (dcdfd = open(dcdname, O_WRONLY|O_LARGEFILE,
00702                                 S_IRUSR|S_IWUSR|S_IRGRP|S_IROTH)) < 0)
00703 #endif
00704         {
00705                 if ( errno != EINTR ) return(DCD_OPENFAILED);
00706         }
00707 
00708         return(dcdfd);
00709 #else
00710         //In the case of multiple output files, each slave has its own output file and
00711         //it needs to be created. --Chao Mei
00712         return open_dcd_write(dcdname);
00713 #endif
00714 }
00715 
00716 /************************************************************************/
00717 /*                                                                      */
00718 /*                              FUNCTION write_dcdstep                  */
00719 /*                                                                      */
00720 /*   INPUTS:                                                            */
00721 /*      fd - file descriptor for the DCD file to write to               */
00722 /*      N - Number of atoms                                             */
00723 /*      X - X coordinates                                               */
00724 /*      Y - Y coordinates                                               */
00725 /*      Z - Z coordinates                                               */
00726 /*  unitcell - a, b, c, alpha, beta, gamma of unit cell */
00727 /*                                                                      */
00728 /*   OUTPUTS:                                                           */
00729 /*      none                                                            */
00730 /*                                                                      */
00731 /*      write_dcdstep writes the coordinates out for a given timestep   */
00732 /*   to the specified DCD file.                                         */
00733 /*                                                                      */
00734 /************************************************************************/
00735 
00736 int write_dcdstep(int fd, int N, float *X, float *Y, float *Z, double *cell)
00737 
00738 {
00739         int32 NSAVC,NSTEP,NFILE;
00740         int32 out_integer;
00741 
00742   /* Unit cell */
00743   if (cell) {
00744     out_integer = 48;
00745     NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00746     NAMD_write(fd, (char *) cell, out_integer);
00747     NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00748   }
00749 
00750   /* Coordinates */
00751   // Note: the value of out_integer wraps for N >= 2^30.
00752         out_integer = N*4;
00753   // Use a separate byte count stored without overflow.
00754   size_t nbytes = ((size_t) N) * 4;
00755 
00756         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00757         NAMD_write(fd, (char *) X, nbytes);
00758         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00759         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00760         NAMD_write(fd, (char *) Y, nbytes);
00761         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00762         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00763         NAMD_write(fd, (char *) Z, nbytes);
00764         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00765 
00766         /* don't update header until after write succeeds */
00767         OFF_T end = LSEEK(fd,0,SEEK_CUR);
00768         LSEEK(fd,NSAVC_POS,SEEK_SET);
00769         READ(fd,(void*) &NSAVC,sizeof(int32));
00770         LSEEK(fd,NSTEP_POS,SEEK_SET);
00771         READ(fd,(void*) &NSTEP,sizeof(int32));
00772         LSEEK(fd,NFILE_POS,SEEK_SET);
00773         READ(fd,(void*) &NFILE,sizeof(int32));
00774         NSTEP += NSAVC;
00775         NFILE += 1;
00776         LSEEK(fd,NSTEP_POS,SEEK_SET);
00777         NAMD_write(fd,(char*) &NSTEP,sizeof(int32));
00778         LSEEK(fd,NFILE_POS,SEEK_SET);
00779         NAMD_write(fd,(char*) &NFILE,sizeof(int32));
00780         LSEEK(fd,end,SEEK_SET);
00781 
00782         return(0);
00783 }
00784 
00785 int write_dcdstep_par_cell(int fd, double *cell){
00786         if (cell) {
00787           int32 out_integer = 48;
00788           NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00789           NAMD_write(fd, (char *) cell, out_integer);
00790           NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00791         }
00792         return 0;
00793 }
00794 
00795 /*   No useful file format description of the DCD format can be found.  */
00796 /*   The closest approximation purporting to be a format description is */
00797 /*   a block of fortran formated statements.  Ultra lame.               */
00798 
00799 /*   Therefore I simply reverse engineered the sequential output.       */
00800 /*    -EJB */
00801 
00802 
00803 
00804 /*   Notionally, one creates the file on the master with header and cell*/
00805 /*   Master then broadcasts a command for the slaves to write.  They    */
00806 /*   contribute to a reduction which roots at the master.  Which then   */
00807 /*   updates the header.                                                */  
00808 /*   Update the X/Y/Z headers for each X,Y and Z fields starting the current */
00809 /*   position of the file                                                    */
00810 int write_dcdstep_par_XYZUnits(int fd, int N)
00811 {
00812   int32 out_integer;
00813 
00814   // number of elements
00815   // Note: the value of out_integer wraps for N >= 2^30.
00816   out_integer = N*sizeof(float);
00817 
00818   // For byte seeking, use a properly sized variable.
00819   OFF_T nbytes = ((OFF_T) N) * sizeof(float);
00820 
00821   NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00822   // seek to the end of each x y z block and write out the count
00823   LSEEK(fd, nbytes, SEEK_CUR);
00824   NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00825 
00826   NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00827   LSEEK(fd, nbytes, SEEK_CUR);
00828   NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00829 
00830   NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00831   LSEEK(fd, nbytes, SEEK_CUR);
00832   NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00833 
00834   return(0);
00835 }
00836 
00837 /*  A function to update the header.                                      */
00838 /*  Once all slaves have written the master updates the header.           */
00839 int update_dcdstep_par_header(int fd)
00840 {
00841         int32 NSAVC,NSTEP,NFILE;
00842         /* don't update header until after write succeeds */
00843         OFF_T end = LSEEK(fd,0,SEEK_CUR);
00844         LSEEK(fd,NSAVC_POS,SEEK_SET);
00845         READ(fd,(void*) &NSAVC,sizeof(int32));
00846         LSEEK(fd,NSTEP_POS,SEEK_SET);
00847         READ(fd,(void*) &NSTEP,sizeof(int32));
00848         LSEEK(fd,NFILE_POS,SEEK_SET);
00849         READ(fd,(void*) &NFILE,sizeof(int32));
00850         NSTEP += NSAVC;
00851         NFILE += 1;
00852         LSEEK(fd,NSTEP_POS,SEEK_SET);
00853         NAMD_write(fd,(char*) &NSTEP,sizeof(int32));
00854         LSEEK(fd,NFILE_POS,SEEK_SET);
00855         NAMD_write(fd,(char*) &NFILE,sizeof(int32));
00856         LSEEK(fd,end,SEEK_SET);
00857 
00858         return(0);
00859 }
00860 
00861 /* Each slave writer has a sequential block of atoms to output */
00862 /* The stream file position needs to be put at the beginning of each */
00863 /* timestep X/Y/Z output                                             */
00864 int write_dcdstep_par_slave(int fd, int parL, int parU, int N, float *X, float *Y, float *Z){
00865 
00866   /* Coordinates for the N elements handled by this writer */
00867         int parN = parU-parL+1;
00868   size_t nbytes = ((size_t)parN) * 4;
00869 
00870         /* x's 1st number of Xs */
00871         /* skip field for the bytes in X, and the first parL atoms in X*/
00872         OFF_T xoffset = sizeof(int)+sizeof(float)*((OFF_T)parL);
00873         LSEEK(fd, xoffset, SEEK_CUR);
00874         NAMD_write(fd, (char *) X, nbytes);
00875 
00876         /* skip field for the bytes in X and Y; */
00877         /* skip the remaining atoms in X at number of (N-1)-(parU+1)+1
00878          * where N-1 is the last atom id, paru+1 is the next atom id. */
00879         /* skip the first parL atoms in Y; */
00880         OFF_T yoffset = 2*sizeof(int)+sizeof(float)*((OFF_T)(N-parU+parL-1));
00881         LSEEK(fd, yoffset, SEEK_CUR);
00882         NAMD_write(fd, (char *) Y, nbytes);
00883 
00884         OFF_T zoffset = yoffset;
00885         LSEEK(fd, zoffset, SEEK_CUR);
00886         NAMD_write(fd, (char *) Z, nbytes);
00887         return(0);
00888 }
00889 
00890 /*****************************************************************************/
00891 /*                                                                           */
00892 /*                              FUNCTION write_dcdheader                     */
00893 /*                                                                           */
00894 /*   INPUTS:                                                                 */
00895 /*      fd - file descriptor for the dcd file                                */
00896 /*      filename - filename for output                                       */
00897 /*      N - Number of atoms                                                  */
00898 /*      NFILE - Number of sets of coordinates                                */
00899 /*      NPRIV - Starting timestep of DCD file - NOT ZERO                     */
00900 /*      NSAVC - Timesteps between DCD saves                                  */
00901 /*      NSTEP - Number of timesteps                                          */
00902 /*      DELTA - length of a timestep                                         */
00903 /*                                                                           */
00904 /*   OUTPUTS:                                                                */
00905 /*      none                                                                 */
00906 /*                                                                           */
00907 /*      This function prints the "header" information to the DCD file.  Since*/
00908 /*   this is duplicating an unformatted binary output from FORTRAN, its ugly.*/
00909 /*   So if you're squimish, don't look.                                      */
00910 /* Whenever you are updating this function by removing or adding fields, please */
00911 /* also update get_dcdheader_size function to reflect the changes! */
00912 /*                                                                           */
00913 /*****************************************************************************/
00914 
00915 int write_dcdheader(int fd, const char *filename, int N, int NFILE, int NPRIV, 
00916                    int NSAVC, int NSTEP, double DELTA, int with_unitcell)
00917 {
00918         int32   out_integer;
00919         float   out_float;
00920         char    title_string[200];
00921         int     user_id;
00922 #ifndef WIN32
00923         struct  passwd *pwbuf;
00924 #endif
00925         time_t  cur_time;
00926         struct  tm *tmbuf;
00927         char    time_str[11];
00928 
00929         out_integer = 84;
00930         NAMD_write(fd, (char *) & out_integer, sizeof(int32));
00931         strcpy(title_string, "CORD");
00932         NAMD_write(fd, title_string, 4);
00933         out_integer = NFILE;  /* located at fpos 8 */
00934         out_integer = 0;  /* ignore the lies */
00935         NAMD_write(fd, (char *) & out_integer, sizeof(int32));
00936         out_integer = NPRIV;  /* located at fpos 12 */
00937         NAMD_write(fd, (char *) & out_integer, sizeof(int32));
00938         out_integer = NSAVC;  /* located at fpos 16 */
00939         NAMD_write(fd, (char *) & out_integer, sizeof(int32));
00940         out_integer = NSTEP;  /* located at fpos 20 */
00941         out_integer = NPRIV - NSAVC;  /* ignore the lies */
00942         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00943         out_integer=0;
00944         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00945         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00946         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00947         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00948         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00949         out_float = DELTA;
00950         NAMD_write(fd, (char *) &out_float, sizeof(float));
00951   out_integer = with_unitcell ? 1 : 0;
00952         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00953   out_integer = 0;
00954         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00955         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00956         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00957         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00958         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00959         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00960         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00961         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00962         out_integer = 24;  // PRETEND TO BE CHARMM24 -JCP
00963         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00964         out_integer = 84;
00965         NAMD_write(fd, (char *) & out_integer, sizeof(int32));
00966 
00967         out_integer = 164;
00968         NAMD_write(fd, (char *) & out_integer, sizeof(int32));
00969         out_integer = 2;
00970         NAMD_write(fd, (char *) & out_integer, sizeof(int32));
00971 
00972         sprintf(title_string, "REMARKS FILENAME=%s CREATED BY NAMD", filename);
00973         pad(title_string, 80);
00974         NAMD_write(fd, title_string, 80);
00975 
00976         char username[100];
00977 #if defined(WIN32) || defined(NO_GETPWUID)
00978         sprintf(username,"Win32");
00979 #else
00980         user_id= (int) getuid();
00981         pwbuf=getpwuid(user_id);
00982         if ( pwbuf ) sprintf(username,"%s",pwbuf->pw_name);
00983         else sprintf(username,"%d",user_id);
00984 #endif
00985         cur_time=time(NULL);
00986         tmbuf=localtime(&cur_time);
00987         strftime(time_str, 10, "%m/%d/%y", tmbuf);
00988 
00989         sprintf(title_string, "REMARKS DATE: %s CREATED BY USER: %s",
00990            time_str, username);
00991         pad(title_string, 80);
00992         NAMD_write(fd, title_string, 80);
00993         out_integer = 164;
00994         NAMD_write(fd, (char *) & out_integer, sizeof(int32));
00995         out_integer = 4;
00996         NAMD_write(fd, (char *) & out_integer, sizeof(int32));
00997         out_integer = N;
00998         NAMD_write(fd, (char *) & out_integer, sizeof(int32));
00999         out_integer = 4;
01000         NAMD_write(fd, (char *) & out_integer, sizeof(int32));
01001 
01002         return(0);
01003 }
01004 
01005 int get_dcdheader_size(){
01006         int headersize = 0;
01007         int totalInt32s = 27; /* 27 writes from out_integer */
01008         int totalChars = 164; /* 3 writes from title_string */
01009         int totalFloats = 1; /* 1 write from out_float */
01010         headersize = sizeof(int32)*totalInt32s+totalChars+sizeof(float)*totalFloats;
01011         return headersize;
01012 }
01013 
01014 /****************************************************************/
01015 /*                                                              */
01016 /*                      FUNCTION close_dcd_read                 */
01017 /*                                                              */
01018 /*   INPUTS:                                                    */
01019 /*      fd - file descriptor to close                           */
01020 /*      num_fixed - the number of fixed atoms                   */
01021 /*      indexes - Array of indexes to be deallocated            */
01022 /*                                                              */
01023 /*   OUTPUTS:                                                   */
01024 /*      the file pointed to by fd is closed and the memory      */
01025 /*   pointed to by indexes is freed if any was allocated        */
01026 /*                                                              */
01027 /*      close_dcd_read closes a dcd file that was opened for    */
01028 /*   reading.  It also deallocated memory used for the indexes  */
01029 /*   used for the free atom list, if there were fixed atoms.    */
01030 /*                                                              */
01031 /****************************************************************/
01032 
01033 void close_dcd_read(int fd, int num_fixed, int *indexes)
01034 
01035 {
01036 #ifdef WIN32
01037         _close(fd);
01038 #else
01039         close(fd);
01040 #endif
01041 
01042         if (num_fixed)
01043         {
01044                 delete [] indexes;
01045         }
01046 }
01047 
01048 /****************************************************************/
01049 /*                                                              */
01050 /*                      FUNCTION close_dcd_write                */
01051 /*                                                              */
01052 /*   INPUTS:                                                    */
01053 /*      fd - file descriptor to close                           */
01054 /*                                                              */
01055 /*   OUTPUTS:                                                   */
01056 /*      the file pointed to by fd                               */
01057 /*                                                              */
01058 /*      close_dcd_write close a dcd file that was opened for    */
01059 /*   writing                                                    */
01060 /*                                                              */
01061 /****************************************************************/
01062 
01063 void close_dcd_write(int fd)
01064 
01065 {
01066 #ifdef WIN32
01067   if ( _close(fd) )
01068 #else
01069   if ( fsync(fd) || close(fd) )
01070 #endif
01071   {
01072     NAMD_err("Error closing DCD file");
01073   }
01074 }
01075 

Generated on Wed Jan 17 01:17:12 2018 for NAMD by  doxygen 1.4.7