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         out_integer = N*4;
00752         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00753         NAMD_write(fd, (char *) X, out_integer);
00754         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00755         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00756         NAMD_write(fd, (char *) Y, out_integer);
00757         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00758         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00759         NAMD_write(fd, (char *) Z, out_integer);
00760         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00761 
00762         /* don't update header until after write succeeds */
00763         OFF_T end = LSEEK(fd,0,SEEK_CUR);
00764         LSEEK(fd,NSAVC_POS,SEEK_SET);
00765         READ(fd,(void*) &NSAVC,sizeof(int32));
00766         LSEEK(fd,NSTEP_POS,SEEK_SET);
00767         READ(fd,(void*) &NSTEP,sizeof(int32));
00768         LSEEK(fd,NFILE_POS,SEEK_SET);
00769         READ(fd,(void*) &NFILE,sizeof(int32));
00770         NSTEP += NSAVC;
00771         NFILE += 1;
00772         LSEEK(fd,NSTEP_POS,SEEK_SET);
00773         NAMD_write(fd,(char*) &NSTEP,sizeof(int32));
00774         LSEEK(fd,NFILE_POS,SEEK_SET);
00775         NAMD_write(fd,(char*) &NFILE,sizeof(int32));
00776         LSEEK(fd,end,SEEK_SET);
00777 
00778         return(0);
00779 }
00780 
00781 int write_dcdstep_par_cell(int fd, double *cell){
00782         if (cell) {
00783           int32 out_integer = 48;
00784           NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00785           NAMD_write(fd, (char *) cell, out_integer);
00786           NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00787         }
00788         return 0;
00789 }
00790 
00791 /*   No useful file format description of the DCD format can be found.  */
00792 /*   The closest approximation purporting to be a format description is */
00793 /*   a block of fortran formated statements.  Ultra lame.               */
00794 
00795 /*   Therefore I simply reverse engineered the sequential output.       */
00796 /*    -EJB */
00797 
00798 
00799 
00800 /*   Notionally, one creates the file on the master with header and cell*/
00801 /*   Master then broadcasts a command for the slaves to write.  They    */
00802 /*   contribute to a reduction which roots at the master.  Which then   */
00803 /*   updates the header.                                                */  
00804 /*   Update the X/Y/Z headers for each X,Y and Z fields starting the current */
00805 /*   position of the file                                                    */
00806 int write_dcdstep_par_XYZUnits(int fd, int N)
00807 {
00808   int32 out_integer;
00809 
00810   // number of elements
00811   out_integer = N*sizeof(float);
00812   NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00813   // seek to the end of each x y z block and write out the count
00814   LSEEK(fd, out_integer, SEEK_CUR);
00815   NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00816 
00817   NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00818   LSEEK(fd, out_integer, SEEK_CUR);
00819   NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00820 
00821   NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00822   LSEEK(fd, out_integer, SEEK_CUR);
00823   NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00824 
00825   return(0);
00826 }
00827 
00828 /*  A function to update the header.                                      */
00829 /*  Once all slaves have written the master updates the header.           */
00830 int update_dcdstep_par_header(int fd)
00831 {
00832         int32 NSAVC,NSTEP,NFILE;
00833         /* don't update header until after write succeeds */
00834         OFF_T end = LSEEK(fd,0,SEEK_CUR);
00835         LSEEK(fd,NSAVC_POS,SEEK_SET);
00836         READ(fd,(void*) &NSAVC,sizeof(int32));
00837         LSEEK(fd,NSTEP_POS,SEEK_SET);
00838         READ(fd,(void*) &NSTEP,sizeof(int32));
00839         LSEEK(fd,NFILE_POS,SEEK_SET);
00840         READ(fd,(void*) &NFILE,sizeof(int32));
00841         NSTEP += NSAVC;
00842         NFILE += 1;
00843         LSEEK(fd,NSTEP_POS,SEEK_SET);
00844         NAMD_write(fd,(char*) &NSTEP,sizeof(int32));
00845         LSEEK(fd,NFILE_POS,SEEK_SET);
00846         NAMD_write(fd,(char*) &NFILE,sizeof(int32));
00847         LSEEK(fd,end,SEEK_SET);
00848 
00849         return(0);
00850 }
00851 
00852 /* Each slave writer has a sequential block of atoms to output */
00853 /* The stream file position needs to be put at the beginning of each */
00854 /* timestep X/Y/Z output                                             */
00855 int write_dcdstep_par_slave(int fd, int parL, int parU, int N, float *X, float *Y, float *Z){
00856 
00857         int parN = parU-parL+1;
00858         int32 out_integer;
00859   /* Coordinates for the N elements handled by this writer */
00860         out_integer = parN*4;
00861 
00862         /* x's 1st number of Xs */
00863         /* skip field for the bytes in X, and the first parL atoms in X*/
00864         OFF_T xoffset = sizeof(int)+sizeof(float)*((OFF_T)parL);
00865         LSEEK(fd, xoffset, SEEK_CUR);
00866         NAMD_write(fd, (char *) X, out_integer);
00867 
00868         /* skip field for the bytes in X and Y; */
00869         /* skip the remaining atoms in X at number of (N-1)-(parU+1)+1
00870          * where N-1 is the last atom id, paru+1 is the next atom id. */
00871         /* skip the first parL atoms in Y; */
00872         OFF_T yoffset = 2*sizeof(int)+sizeof(float)*((OFF_T)(N-parU+parL-1));
00873         LSEEK(fd, yoffset, SEEK_CUR);
00874         NAMD_write(fd, (char *) Y, out_integer);
00875 
00876         OFF_T zoffset = yoffset;
00877         LSEEK(fd, zoffset, SEEK_CUR);
00878         NAMD_write(fd, (char *) Z, out_integer);
00879         return(0);
00880 }
00881 
00882 /*****************************************************************************/
00883 /*                                                                           */
00884 /*                              FUNCTION write_dcdheader                     */
00885 /*                                                                           */
00886 /*   INPUTS:                                                                 */
00887 /*      fd - file descriptor for the dcd file                                */
00888 /*      filename - filename for output                                       */
00889 /*      N - Number of atoms                                                  */
00890 /*      NFILE - Number of sets of coordinates                                */
00891 /*      NPRIV - Starting timestep of DCD file - NOT ZERO                     */
00892 /*      NSAVC - Timesteps between DCD saves                                  */
00893 /*      NSTEP - Number of timesteps                                          */
00894 /*      DELTA - length of a timestep                                         */
00895 /*                                                                           */
00896 /*   OUTPUTS:                                                                */
00897 /*      none                                                                 */
00898 /*                                                                           */
00899 /*      This function prints the "header" information to the DCD file.  Since*/
00900 /*   this is duplicating an unformatted binary output from FORTRAN, its ugly.*/
00901 /*   So if you're squimish, don't look.                                      */
00902 /* Whenever you are updating this function by removing or adding fields, please */
00903 /* also update get_dcdheader_size function to reflect the changes! */
00904 /*                                                                           */
00905 /*****************************************************************************/
00906 
00907 int write_dcdheader(int fd, const char *filename, int N, int NFILE, int NPRIV, 
00908                    int NSAVC, int NSTEP, double DELTA, int with_unitcell)
00909 {
00910         int32   out_integer;
00911         float   out_float;
00912         char    title_string[200];
00913         int     user_id;
00914 #ifndef WIN32
00915         struct  passwd *pwbuf;
00916 #endif
00917         time_t  cur_time;
00918         struct  tm *tmbuf;
00919         char    time_str[11];
00920 
00921         out_integer = 84;
00922         NAMD_write(fd, (char *) & out_integer, sizeof(int32));
00923         strcpy(title_string, "CORD");
00924         NAMD_write(fd, title_string, 4);
00925         out_integer = NFILE;  /* located at fpos 8 */
00926         out_integer = 0;  /* ignore the lies */
00927         NAMD_write(fd, (char *) & out_integer, sizeof(int32));
00928         out_integer = NPRIV;  /* located at fpos 12 */
00929         NAMD_write(fd, (char *) & out_integer, sizeof(int32));
00930         out_integer = NSAVC;  /* located at fpos 16 */
00931         NAMD_write(fd, (char *) & out_integer, sizeof(int32));
00932         out_integer = NSTEP;  /* located at fpos 20 */
00933         out_integer = NPRIV - NSAVC;  /* ignore the lies */
00934         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00935         out_integer=0;
00936         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00937         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00938         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00939         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00940         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00941         out_float = DELTA;
00942         NAMD_write(fd, (char *) &out_float, sizeof(float));
00943   out_integer = with_unitcell ? 1 : 0;
00944         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00945   out_integer = 0;
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         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00950         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00951         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00952         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00953         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00954         out_integer = 24;  // PRETEND TO BE CHARMM24 -JCP
00955         NAMD_write(fd, (char *) &out_integer, sizeof(int32));
00956         out_integer = 84;
00957         NAMD_write(fd, (char *) & out_integer, sizeof(int32));
00958 
00959         out_integer = 164;
00960         NAMD_write(fd, (char *) & out_integer, sizeof(int32));
00961         out_integer = 2;
00962         NAMD_write(fd, (char *) & out_integer, sizeof(int32));
00963 
00964         sprintf(title_string, "REMARKS FILENAME=%s CREATED BY NAMD", filename);
00965         pad(title_string, 80);
00966         NAMD_write(fd, title_string, 80);
00967 
00968         char username[100];
00969 #if defined(WIN32) || defined(NO_GETPWUID)
00970         sprintf(username,"Win32");
00971 #else
00972         user_id= (int) getuid();
00973         pwbuf=getpwuid(user_id);
00974         if ( pwbuf ) sprintf(username,"%s",pwbuf->pw_name);
00975         else sprintf(username,"%d",user_id);
00976 #endif
00977         cur_time=time(NULL);
00978         tmbuf=localtime(&cur_time);
00979         strftime(time_str, 10, "%m/%d/%y", tmbuf);
00980 
00981         sprintf(title_string, "REMARKS DATE: %s CREATED BY USER: %s",
00982            time_str, username);
00983         pad(title_string, 80);
00984         NAMD_write(fd, title_string, 80);
00985         out_integer = 164;
00986         NAMD_write(fd, (char *) & out_integer, sizeof(int32));
00987         out_integer = 4;
00988         NAMD_write(fd, (char *) & out_integer, sizeof(int32));
00989         out_integer = N;
00990         NAMD_write(fd, (char *) & out_integer, sizeof(int32));
00991         out_integer = 4;
00992         NAMD_write(fd, (char *) & out_integer, sizeof(int32));
00993 
00994         return(0);
00995 }
00996 
00997 int get_dcdheader_size(){
00998         int headersize = 0;
00999         int totalInt32s = 27; /* 27 writes from out_integer */
01000         int totalChars = 164; /* 3 writes from title_string */
01001         int totalFloats = 1; /* 1 write from out_float */
01002         headersize = sizeof(int32)*totalInt32s+totalChars+sizeof(float)*totalFloats;
01003         return headersize;
01004 }
01005 
01006 /****************************************************************/
01007 /*                                                              */
01008 /*                      FUNCTION close_dcd_read                 */
01009 /*                                                              */
01010 /*   INPUTS:                                                    */
01011 /*      fd - file descriptor to close                           */
01012 /*      num_fixed - the number of fixed atoms                   */
01013 /*      indexes - Array of indexes to be deallocated            */
01014 /*                                                              */
01015 /*   OUTPUTS:                                                   */
01016 /*      the file pointed to by fd is closed and the memory      */
01017 /*   pointed to by indexes is freed if any was allocated        */
01018 /*                                                              */
01019 /*      close_dcd_read closes a dcd file that was opened for    */
01020 /*   reading.  It also deallocated memory used for the indexes  */
01021 /*   used for the free atom list, if there were fixed atoms.    */
01022 /*                                                              */
01023 /****************************************************************/
01024 
01025 void close_dcd_read(int fd, int num_fixed, int *indexes)
01026 
01027 {
01028 #ifdef WIN32
01029         _close(fd);
01030 #else
01031         close(fd);
01032 #endif
01033 
01034         if (num_fixed)
01035         {
01036                 delete [] indexes;
01037         }
01038 }
01039 
01040 /****************************************************************/
01041 /*                                                              */
01042 /*                      FUNCTION close_dcd_write                */
01043 /*                                                              */
01044 /*   INPUTS:                                                    */
01045 /*      fd - file descriptor to close                           */
01046 /*                                                              */
01047 /*   OUTPUTS:                                                   */
01048 /*      the file pointed to by fd                               */
01049 /*                                                              */
01050 /*      close_dcd_write close a dcd file that was opened for    */
01051 /*   writing                                                    */
01052 /*                                                              */
01053 /****************************************************************/
01054 
01055 void close_dcd_write(int fd)
01056 
01057 {
01058 #ifdef WIN32
01059   if ( _close(fd) )
01060 #else
01061   if ( fsync(fd) || close(fd) )
01062 #endif
01063   {
01064     NAMD_err("Error closing DCD file");
01065   }
01066 }
01067 

Generated on Sat Nov 18 01:17:13 2017 for NAMD by  doxygen 1.4.7