NAMD
dcdlib.C
Go to the documentation of this file.
1 
7 /*
8  dcdlib contains C routines for reading and writing binary DCD
9  files. The output format of these files is based on binary FORTRAN
10  output, so its pretty ugly. If you are squimish, don't look!
11 */
12 
13 #include "dcdlib.h"
14 
15 #ifndef OUTPUT_SINGLE_FILE
16 #define OUTPUT_SINGLE_FILE 1
17 // including Output.h causes compile error on BlueGeneQ
18 #endif
19 
20 #ifdef WIN32
21 #define access(PATH,MODE) _access(PATH,00)
22 #endif
23 
24 #define NAMD_write NAMD_write64
25 // same as write, only does error checking internally
26 void NAMD_write(int fd, const char *buf, size_t count) {
27  while ( count ) {
28 #if defined(WIN32) && !defined(__CYGWIN__)
29  long retval = _write(fd,buf,count);
30 #else
31  ssize_t retval = write(fd,buf,count);
32 #endif
33  if ( retval < 0 && errno == EINTR ) retval = 0;
34  if ( retval < 0 ) NAMD_die(strerror(errno));
35  if ( retval > count ) NAMD_bug("extra bytes written in NAMD_write64()");
36  buf += retval;
37  count -= retval;
38  }
39 }
40 
41 #ifdef WIN32
42 #define LSEEK _lseeki64
43 #define READ _read
44 #else
45 #define LSEEK lseek
46 #define READ read
47 #endif
48 
49 OFF_T NAMD_seek(int file, OFF_T offset, int whence) {
50  OFF_T retval = LSEEK(file, offset, whence);
51  if ( retval < 0 ) NAMD_err("seek failed while writing DCD file");
52  if ( whence == SEEK_SET && retval != offset ) {
53  char buf[256];
54  sprintf(buf, "seek failed while writing DCD file: SEEK_SET %lld returned %lld\n", offset, retval);
55  NAMD_die(buf);
56  }
57  return retval;
58 }
59 
60 #undef LSEEK
61 #define LSEEK NAMD_seek
62 
63 #ifndef O_LARGEFILE
64 #define O_LARGEFILE 0x0
65 #endif
66 
67 
68 /************************************************************************/
69 /* */
70 /* MACRO CHECK_FREAD */
71 /* */
72 /* CHECK_FREAD is used to check the status of a file read. It */
73 /* is passed the return code from read and a string to print out if */
74 /* an error is detected. If an error is found, an error message is */
75 /* printed out and the program terminates. This was made into a */
76 /* macro because it had to be done over and over and over . . . */
77 /* */
78 /************************************************************************/
79 
80 #define CHECK_FREAD(X, msg) if (X==-1) \
81  { \
82  return(DCD_BADREAD); \
83  }
84 
85 /************************************************************************/
86 /* */
87 /* MACRO CHECK_FEOF */
88 /* */
89 /* CHECK_FEOF checks for an abnormal EOF on a file read. It is */
90 /* passed the return status from read and a message to print on error.*/
91 /* if an EOF is detected, the error message is printed and the program*/
92 /* terminates. This is used for reads that should find something in */
93 /* the file. */
94 /* */
95 /************************************************************************/
96 
97 #define CHECK_FEOF(X, msg) if (X==0) \
98  { \
99  return(DCD_BADEOF); \
100  }
101 
102 
103 /* GLOBAL VARIABLES */
104 
105 void pad(char *s, int len)
106 {
107  int curlen;
108  int i;
109 
110  curlen=strlen(s);
111 
112  if (curlen>len)
113  {
114  s[len]='\0';
115  return;
116  }
117 
118  for (i=curlen; i<len; i++)
119  {
120  s[i]=' ';
121  }
122 
123  s[i]='\0';
124 }
125 
126 /************************************************************************/
127 /* */
128 /* FUNCTION open_dcd_read */
129 /* */
130 /* INPUTS: */
131 /* filename - filename to read in */
132 /* */
133 /* OUTPUTS: */
134 /* an open filedesriptor (integer) is returned if the call is */
135 /* successful, otherwise a negative number is returned */
136 /* */
137 /* open_dcd_read opens a dcd file for input. It first does a check*/
138 /* to see if the file really exists. If the file does not exist, */
139 /* a DCD_DNE is returned, if the file exists but can' be opened for */
140 /* some reason, a DCD_OPENFAILED is returned. If the open is */
141 /* successful, the file descriptor is returned. */
142 /* */
143 /************************************************************************/
144 
145 int open_dcd_read(char *filename)
146 
147 {
148  int dcdfd; /* file descriptor for dcd file */
149 
150  /* Do a stat just to see if the file really exists */
151  if (access(filename, F_OK) != 0)
152  {
153  if (errno == ENOENT)
154  {
155  return(DCD_DNE);
156  }
157  }
158 
159  /* Try and open the file */
160 #ifdef WIN32
161  dcdfd=_open(filename, O_RDONLY|O_BINARY|O_LARGEFILE);
162 #else
163  dcdfd=open(filename, O_RDONLY|O_LARGEFILE);
164 #endif
165 
166  if (dcdfd == -1)
167  {
168  return(DCD_OPENFAILED);
169  }
170 
171  return(dcdfd);
172 }
173 
174 /****************************************************************/
175 /* */
176 /* FUNCTION read_dcdheader */
177 /* */
178 /* INPUTS: */
179 /* fd - file descriptor for the dcd file */
180 /* N - Number of atoms */
181 /* NSET - Number of sets of coordinates */
182 /* ISTART - Starting timestep of DCD file */
183 /* NSAVC - Timesteps between DCD saves */
184 /* DELTA - length of a timestep */
185 /* */
186 /* OUTPUTS: */
187 /* N, NSET, ISTART, NSAVC, and DELTA are returned populated*/
188 /* a 0 is returned if the call is successful, or a negative */
189 /* number if errors are detected */
190 /* */
191 /* read_header reads in the header from a DCD file and */
192 /* returns the timestep size and the number of atoms in the */
193 /* system. A 0 is returned if the header is successfully */
194 /* read, or a negative error code is returned otherwise. */
195 /* */
196 /****************************************************************/
197 
198 #if 0
199 int read_dcdheader(int fd, int *N, int *NSET, int *ISTART,
200  int *NSAVC, double *DELTA, int *NAMNF,
201  int **FREEINDEXES)
202 
203 {
204  int32 input_integer; /* Integer buffer space */
205  float input_float;
206  char bigbuf[256]; /* A large string buffer */
207  int ret_val; /* Return value from read */
208  int i; /* Loop counter */
209  char HDR[5]; /* Header = "CORD" */
210  int32 I;
211  int32 NTITLE;
212 
213  /* First thing in the file should be an 84 */
214  ret_val = read(fd, &input_integer, sizeof(int32));
215 
216  CHECK_FREAD(ret_val, "reading first int from dcd file");
217  CHECK_FEOF(ret_val, "reading first int from dcd file");
218 
219  if (input_integer != 84)
220  {
221  return(DCD_BADFORMAT);
222  }
223 
224  /* Read in the string "CORD" */
225  ret_val = read(fd, HDR, 4);
226 
227  CHECK_FREAD(ret_val, "reading CORD from dcd file");
228  CHECK_FEOF(ret_val, "reading CORD from dcd file");
229 
230  HDR[ret_val]='\0';
231 
232  if (strcmp(HDR, "CORD") != 0)
233  {
234  return(DCD_BADFORMAT);
235  }
236 
237  /* Read in the number of Sets of coordinates, NSET */
238  ret_val = read(fd, &input_integer, sizeof(int32));
239  *NSET = input_integer;
240 
241  CHECK_FREAD(ret_val, "reading NSET");
242  CHECK_FEOF(ret_val, "reading NSET");
243 
244  /* Read in ISTART, the starting timestep */
245  ret_val = read(fd, &input_integer, sizeof(int32));
246  *ISTART = input_integer;
247 
248  CHECK_FREAD(ret_val, "reading ISTART");
249  CHECK_FEOF(ret_val, "reading ISTART");
250 
251  /* Read in NSAVC, the number of timesteps between */
252  /* dcd saves */
253  ret_val = read(fd, &input_integer, sizeof(int32));
254  *NSAVC = input_integer;
255 
256  CHECK_FREAD(ret_val, "reading NSAVC");
257  CHECK_FEOF(ret_val, "reading NSAVC");
258 
259 
260  /* Skip blank integers */
261  for (i=0; i<5; i++)
262  {
263  ret_val = read(fd, &input_integer, sizeof(int32));
264 
265  CHECK_FREAD(ret_val, "reading I");
266  CHECK_FEOF(ret_val, "reading I");
267  }
268 
269  /* Read NAMNF, the number of free atoms */
270  ret_val = read(fd, &input_integer, sizeof(int32));
271  *NAMNF = input_integer;
272 
273  CHECK_FREAD(ret_val, "reading NAMNF");
274  CHECK_FEOF(ret_val, "reading NAMNF");
275 
276  /* Read in the timestep, DELTA */
277  ret_val = read(fd, &input_float, sizeof(float));
278 
279  CHECK_FREAD(ret_val, "reading DELTA");
280  CHECK_FEOF(ret_val, "reading DELTA");
281  *DELTA = input_float;
282 
283  /* Skip blank integers */
284  for (i=0; i<10; i++)
285  {
286  ret_val = read(fd, &I, sizeof(int32));
287 
288  CHECK_FREAD(ret_val, "reading I");
289  CHECK_FEOF(ret_val, "reading I");
290  }
291 
292  /* Get the end size of the first block */
293  ret_val = read(fd, &input_integer, sizeof(int32));
294 
295  CHECK_FREAD(ret_val, "reading second 84 from dcd file");
296  CHECK_FEOF(ret_val, "reading second 84 from dcd file");
297 
298  if (input_integer != 84)
299  {
300  return(DCD_BADFORMAT);
301  }
302 
303  /* Read in the size of the next block */
304  ret_val = read(fd, &input_integer, sizeof(int32));
305 
306  CHECK_FREAD(ret_val, "reading size of title block");
307  CHECK_FEOF(ret_val, "reading size of title block");
308 
309  if ( ((input_integer-4)%80) == 0)
310  {
311  /* Read NTITLE, the number of 80 characeter */
312  /* title strings there are */
313  ret_val = read(fd, &NTITLE, sizeof(int32));
314 
315  CHECK_FREAD(ret_val, "reading NTITLE");
316  CHECK_FEOF(ret_val, "reading NTITLE");
317 
318  for (i=0; i<NTITLE; i++)
319  {
320  ret_val = read(fd, bigbuf, 80);
321 
322  CHECK_FREAD(ret_val, "reading TITLE");
323  CHECK_FEOF(ret_val, "reading TITLE");
324  }
325 
326  /* Get the ending size for this block */
327  ret_val = read(fd, &input_integer, sizeof(int32));
328 
329  CHECK_FREAD(ret_val, "reading size of title block");
330  CHECK_FEOF(ret_val, "reading size of title block");
331  }
332  else
333  {
334  return(DCD_BADFORMAT);
335  }
336 
337  /* Read in an 4 */
338  ret_val = read(fd, &input_integer, sizeof(int32));
339 
340  CHECK_FREAD(ret_val, "reading an 4");
341  CHECK_FEOF(ret_val, "reading an 4");
342 
343  if (input_integer != 4)
344  {
345  return(DCD_BADFORMAT);
346  }
347 
348  /* Read in the number of atoms */
349  ret_val = read(fd, &input_integer, sizeof(int32));
350  *N = input_integer;
351 
352  CHECK_FREAD(ret_val, "reading number of atoms");
353  CHECK_FEOF(ret_val, "reading number of atoms");
354 
355  /* Read in an 4 */
356  ret_val = read(fd, &input_integer, sizeof(int32));
357 
358  CHECK_FREAD(ret_val, "reading an 4");
359  CHECK_FEOF(ret_val, "reading an 4");
360 
361  if (input_integer != 4)
362  {
363  return(DCD_BADFORMAT);
364  }
365 
366  if (*NAMNF != 0)
367  {
368  (*FREEINDEXES) = new int[(*N)-(*NAMNF)];
369  int32 *freeindexes32 = new int32[(*N)-(*NAMNF)];
370 
371  if (*FREEINDEXES == NULL || freeindexes32 == NULL)
372  return(DCD_BADMALLOC);
373 
374  /* Read in an size */
375  ret_val = read(fd, &input_integer, sizeof(int32));
376 
377  CHECK_FREAD(ret_val, "reading size of index array");
378  CHECK_FEOF(ret_val, "reading size of index array");
379 
380  if (input_integer != ((*N)-(*NAMNF))*4)
381  {
382  return(DCD_BADFORMAT);
383  }
384 
385  ret_val = read(fd, freeindexes32, ((*N)-(*NAMNF))*sizeof(int32));
386 
387  CHECK_FREAD(ret_val, "reading size of index array");
388  CHECK_FEOF(ret_val, "reading size of index array");
389 
390  for (i=0; i<((*N)-(*NAMNF)); ++i)
391  (*FREEINDEXES)[i] = freeindexes32[i];
392 
393  delete [] freeindexes32;
394 
395  ret_val = read(fd, &input_integer, sizeof(int32));
396 
397  CHECK_FREAD(ret_val, "reading size of index array");
398  CHECK_FEOF(ret_val, "reading size of index array");
399 
400  if (input_integer != ((*N)-(*NAMNF))*4)
401  {
402  return(DCD_BADFORMAT);
403  }
404  }
405 
406  return(0);
407 }
408 
409 /************************************************************************/
410 /* */
411 /* FUNCTION read_dcdstep */
412 /* */
413 /* INPUTS: */
414 /* fd - file descriptor to use */
415 /* N - Number of atoms */
416 /* X - array of X coordinates */
417 /* Y - array of Y coordinates */
418 /* Z - array of Z coordinates */
419 /* num_fixed - Number of fixed atoms */
420 /* first - flag 1->first time called */
421 /* indexes - indexes for free atoms */
422 /* */
423 /* OUTPUTS: */
424 /* read step populates the arrays X, Y, Z and returns a 0 if the */
425 /* read is successful. If the read fails, a negative error code is */
426 /* returned. */
427 /* */
428 /* read_step reads in the coordinates from one step. It places */
429 /* these coordinates into the arrays X, Y, and Z. */
430 /* */
431 /************************************************************************/
432 
433 int read_dcdstep(int fd, int N, float *X, float *Y, float *Z, int num_fixed,
434  int first, int *indexes)
435 
436 {
437  int ret_val; /* Return value from read */
438  int32 input_integer; /* Integer buffer space */
439  int i; /* Loop counter */
440  static float *tmpX;
441 
442  if (first && num_fixed)
443  {
444  tmpX = new float[N-num_fixed];
445 
446  if (tmpX==NULL)
447  {
448  return(DCD_BADMALLOC);
449  }
450  }
451 
452  /* Get the first size from the file */
453  ret_val = read(fd, &input_integer, sizeof(int32));
454 
455  CHECK_FREAD(ret_val, "reading number of atoms at begining of step");
456 
457  /* See if we've reached the end of the file */
458  if (ret_val == 0)
459  {
460  delete [] tmpX;
461 
462  return(-1);
463  }
464 
465  if ( (num_fixed==0) || first)
466  {
467  if (input_integer != 4*N)
468  {
469  return(DCD_BADFORMAT);
470  }
471 
472  ret_val = read(fd, X, 4*N);
473 
474  CHECK_FREAD(ret_val, "reading X array");
475  CHECK_FEOF(ret_val, "reading X array");
476 
477  ret_val = read(fd, &input_integer, sizeof(int32));
478 
479  CHECK_FREAD(ret_val, "reading number of atoms after X array");
480  CHECK_FEOF(ret_val, "reading number of atoms after X array");
481 
482  if (input_integer != 4*N)
483  {
484  return(DCD_BADFORMAT);
485  }
486 
487  ret_val = read(fd, &input_integer, sizeof(int32));
488 
489  CHECK_FREAD(ret_val, "reading number of atoms after X array");
490  CHECK_FEOF(ret_val, "reading number of atoms after X array");
491 
492  if (input_integer != 4*N)
493  {
494  return(DCD_BADFORMAT);
495  }
496 
497  ret_val = read(fd, Y, 4*N);
498 
499  CHECK_FREAD(ret_val, "reading Y array");
500  CHECK_FEOF(ret_val, "reading Y array");
501 
502  ret_val = read(fd, &input_integer, sizeof(int32));
503 
504  CHECK_FREAD(ret_val, "reading number of atoms after Y array");
505  CHECK_FEOF(ret_val, "reading number of atoms after Y array");
506 
507  if (input_integer != 4*N)
508  {
509  return(DCD_BADFORMAT);
510  }
511 
512  ret_val = read(fd, &input_integer, sizeof(int32));
513 
514  CHECK_FREAD(ret_val, "reading number of atoms after Y array");
515  CHECK_FEOF(ret_val, "reading number of atoms after Y array");
516 
517  if (input_integer != 4*N)
518  {
519  return(DCD_BADFORMAT);
520  }
521 
522  ret_val = read(fd, Z, 4*N);
523 
524  CHECK_FREAD(ret_val, "reading Z array");
525  CHECK_FEOF(ret_val, "reading Z array");
526 
527  ret_val = read(fd, &input_integer, sizeof(int32));
528 
529  CHECK_FREAD(ret_val, "reading number of atoms after Z array");
530  CHECK_FEOF(ret_val, "reading number of atoms after Z array");
531 
532  if (input_integer != 4*N)
533  {
534  return(DCD_BADFORMAT);
535  }
536  }
537  else
538  {
539  if (input_integer != 4*(N-num_fixed))
540  {
541  return(DCD_BADFORMAT);
542  }
543 
544  ret_val = read(fd, tmpX, 4*(N-num_fixed));
545 
546  CHECK_FREAD(ret_val, "reading Xtmp array");
547  CHECK_FEOF(ret_val, "reading Xtmp array");
548 
549  for (i=0; i<N-num_fixed; i++)
550  {
551  X[indexes[i]-1]=tmpX[i];
552  }
553 
554  ret_val = read(fd, &input_integer, sizeof(int32));
555 
556  CHECK_FREAD(ret_val, "reading number of atoms after X array");
557  CHECK_FEOF(ret_val, "reading number of atoms after X array");
558 
559  if (input_integer != 4*(N-num_fixed))
560  {
561  return(DCD_BADFORMAT);
562  }
563 
564  ret_val = read(fd, &input_integer, sizeof(int32));
565 
566  if (input_integer != 4*(N-num_fixed))
567  {
568  return(DCD_BADFORMAT);
569  }
570 
571  ret_val = read(fd, tmpX, 4*(N-num_fixed));
572 
573  CHECK_FREAD(ret_val, "reading Xtmp array");
574  CHECK_FEOF(ret_val, "reading Xtmp array");
575 
576  for (i=0; i<N-num_fixed; i++)
577  {
578  Y[indexes[i]-1]=tmpX[i];
579  }
580 
581  ret_val = read(fd, &input_integer, sizeof(int32));
582 
583  CHECK_FREAD(ret_val, "reading number of atoms after Y array");
584  CHECK_FEOF(ret_val, "reading number of atoms after Y array");
585 
586  if (input_integer != 4*(N-num_fixed))
587  {
588  return(DCD_BADFORMAT);
589  }
590 
591  ret_val = read(fd, &input_integer, sizeof(int32));
592 
593  if (input_integer != 4*(N-num_fixed))
594  {
595  return(DCD_BADFORMAT);
596  }
597 
598  ret_val = read(fd, tmpX, 4*(N-num_fixed));
599 
600  CHECK_FREAD(ret_val, "reading Xtmp array");
601  CHECK_FEOF(ret_val, "reading Xtmp array");
602 
603  for (i=0; i<N-num_fixed; i++)
604  {
605  Z[indexes[i]-1]=tmpX[i];
606  }
607 
608  ret_val = read(fd, &input_integer, sizeof(int32));
609 
610  CHECK_FREAD(ret_val, "reading number of atoms after Z array");
611  CHECK_FEOF(ret_val, "reading number of atoms after Z array");
612 
613  if (input_integer != 4*(N-num_fixed))
614  {
615  return(DCD_BADFORMAT);
616  }
617  }
618 
619  return(0);
620 }
621 #endif
622 
623 #ifndef OUTPUT_SINGLE_FILE
624 #error OUTPUT_SINGLE_FILE not defined!
625 #endif
626 
627 #if OUTPUT_SINGLE_FILE
628 #define NFILE_POS ((OFF_T) 8)
629 #define NPRIV_POS ((OFF_T) 12)
630 #define NSAVC_POS ((OFF_T) 16)
631 #define NSTEP_POS ((OFF_T) 20)
632 #else
633 //Need to consider extra fields:
634 //1. magic number (int32)
635 //2. file version (float)
636 //3. number of files that contain portion of data (int32)
637 //So the total offset is 12 bytes.
638 #define NFILE_POS ((OFF_T) 20)
639 #define NPRIV_POS ((OFF_T) 24)
640 #define NSAVC_POS ((OFF_T) 28)
641 #define NSTEP_POS ((OFF_T) 32)
642 #endif
643 
644 /*********************************************************************/
645 /* */
646 /* FUNCTION open_dcd_write */
647 /* */
648 /* INPUTS: */
649 /* dcdfile - Name of the dcd file */
650 /* */
651 /* OUTPUTS: */
652 /* returns an open file descriptor for writing */
653 /* */
654 /* This function will open a dcd file for writing. It takes */
655 /* the filename to open as its only argument. It will return a */
656 /* valid file descriptor if successful or DCD_OPENFAILED if the */
657 /* open fails for some reason. If the file specifed already */
658 /* exists, it is renamed by appending .BAK to it. */
659 /* */
660 /*********************************************************************/
661 
662 int open_dcd_write(const char *dcdname)
663 
664 {
665  int dcdfd;
666  NAMD_backup_file(dcdname,".BAK");
667 
668 #ifdef WIN32
669  while ( (dcdfd = _open(dcdname, O_RDWR|O_CREAT|O_EXCL|O_BINARY|O_LARGEFILE,
670  _S_IREAD|_S_IWRITE)) < 0)
671 #else
672 #ifdef NAMD_NO_O_EXCL
673  while ( (dcdfd = open(dcdname, O_RDWR|O_CREAT|O_TRUNC|O_LARGEFILE,
674 #else
675  while ( (dcdfd = open(dcdname, O_RDWR|O_CREAT|O_EXCL|O_LARGEFILE,
676 #endif
677  S_IRUSR|S_IWUSR|S_IRGRP|S_IROTH)) < 0)
678 #endif
679  {
680  if ( errno != EINTR ) return(DCD_OPENFAILED);
681  }
682 
683  return(dcdfd);
684 }
685 
686 // non master nodes open an existing file
687 // master will have taken care of backup and creation before the
688 // broadcast which triggers the slaves
689 int open_dcd_write_par_slave(char *dcdname)
690 
691 {
692 #if OUTPUT_SINGLE_FILE
693  //In the case of single file, the dcd output by slaves has been created
694  //by the master, so the dcd file doesn't need to be created again and
695  //backed up. --Chao Mei
696  int dcdfd;
697 #ifdef WIN32
698  while ( (dcdfd = _open(dcdname, O_WRONLY|O_BINARY|O_LARGEFILE,
699  _S_IREAD|_S_IWRITE)) < 0)
700 #else
701  while ( (dcdfd = open(dcdname, O_WRONLY|O_LARGEFILE,
702  S_IRUSR|S_IWUSR|S_IRGRP|S_IROTH)) < 0)
703 #endif
704  {
705  if ( errno != EINTR ) return(DCD_OPENFAILED);
706  }
707 
708  return(dcdfd);
709 #else
710  //In the case of multiple output files, each slave has its own output file and
711  //it needs to be created. --Chao Mei
712  return open_dcd_write(dcdname);
713 #endif
714 }
715 
716 /************************************************************************/
717 /* */
718 /* FUNCTION write_dcdstep */
719 /* */
720 /* INPUTS: */
721 /* fd - file descriptor for the DCD file to write to */
722 /* N - Number of atoms */
723 /* X - X coordinates */
724 /* Y - Y coordinates */
725 /* Z - Z coordinates */
726 /* unitcell - a, b, c, alpha, beta, gamma of unit cell */
727 /* */
728 /* OUTPUTS: */
729 /* none */
730 /* */
731 /* write_dcdstep writes the coordinates out for a given timestep */
732 /* to the specified DCD file. */
733 /* */
734 /************************************************************************/
735 
736 int write_dcdstep(int fd, int N, float *X, float *Y, float *Z, double *cell)
737 
738 {
739  int32 NSAVC,NSTEP,NFILE;
740  int32 out_integer;
741 
742  /* Unit cell */
743  if (cell) {
744  out_integer = 48;
745  NAMD_write(fd, (char *) &out_integer, sizeof(int32));
746  NAMD_write(fd, (char *) cell, out_integer);
747  NAMD_write(fd, (char *) &out_integer, sizeof(int32));
748  }
749 
750  /* Coordinates */
751  // Note: the value of out_integer wraps for N >= 2^30.
752  out_integer = N*4;
753  // Use a separate byte count stored without overflow.
754  size_t nbytes = ((size_t) N) * 4;
755 
756  NAMD_write(fd, (char *) &out_integer, sizeof(int32));
757  NAMD_write(fd, (char *) X, nbytes);
758  NAMD_write(fd, (char *) &out_integer, sizeof(int32));
759  NAMD_write(fd, (char *) &out_integer, sizeof(int32));
760  NAMD_write(fd, (char *) Y, nbytes);
761  NAMD_write(fd, (char *) &out_integer, sizeof(int32));
762  NAMD_write(fd, (char *) &out_integer, sizeof(int32));
763  NAMD_write(fd, (char *) Z, nbytes);
764  NAMD_write(fd, (char *) &out_integer, sizeof(int32));
765 
766  /* don't update header until after write succeeds */
767  OFF_T end = LSEEK(fd,0,SEEK_CUR);
768  LSEEK(fd,NSAVC_POS,SEEK_SET);
769  READ(fd,(void*) &NSAVC,sizeof(int32));
770  LSEEK(fd,NSTEP_POS,SEEK_SET);
771  READ(fd,(void*) &NSTEP,sizeof(int32));
772  LSEEK(fd,NFILE_POS,SEEK_SET);
773  READ(fd,(void*) &NFILE,sizeof(int32));
774  NSTEP += NSAVC;
775  NFILE += 1;
776  LSEEK(fd,NSTEP_POS,SEEK_SET);
777  NAMD_write(fd,(char*) &NSTEP,sizeof(int32));
778  LSEEK(fd,NFILE_POS,SEEK_SET);
779  NAMD_write(fd,(char*) &NFILE,sizeof(int32));
780  LSEEK(fd,end,SEEK_SET);
781 
782  return(0);
783 }
784 
785 int write_dcdstep_par_cell(int fd, double *cell){
786  if (cell) {
787  int32 out_integer = 48;
788  NAMD_write(fd, (char *) &out_integer, sizeof(int32));
789  NAMD_write(fd, (char *) cell, out_integer);
790  NAMD_write(fd, (char *) &out_integer, sizeof(int32));
791  }
792  return 0;
793 }
794 
795 /* No useful file format description of the DCD format can be found. */
796 /* The closest approximation purporting to be a format description is */
797 /* a block of fortran formated statements. Ultra lame. */
798 
799 /* Therefore I simply reverse engineered the sequential output. */
800 /* -EJB */
801 
802 
803 
804 /* Notionally, one creates the file on the master with header and cell*/
805 /* Master then broadcasts a command for the slaves to write. They */
806 /* contribute to a reduction which roots at the master. Which then */
807 /* updates the header. */
808 /* Update the X/Y/Z headers for each X,Y and Z fields starting the current */
809 /* position of the file */
810 int write_dcdstep_par_XYZUnits(int fd, int N)
811 {
812  int32 out_integer;
813 
814  // number of elements
815  // Note: the value of out_integer wraps for N >= 2^30.
816  out_integer = N*sizeof(float);
817 
818  // For byte seeking, use a properly sized variable.
819  OFF_T nbytes = ((OFF_T) N) * sizeof(float);
820 
821  NAMD_write(fd, (char *) &out_integer, sizeof(int32));
822  // seek to the end of each x y z block and write out the count
823  LSEEK(fd, nbytes, SEEK_CUR);
824  NAMD_write(fd, (char *) &out_integer, sizeof(int32));
825 
826  NAMD_write(fd, (char *) &out_integer, sizeof(int32));
827  LSEEK(fd, nbytes, SEEK_CUR);
828  NAMD_write(fd, (char *) &out_integer, sizeof(int32));
829 
830  NAMD_write(fd, (char *) &out_integer, sizeof(int32));
831  LSEEK(fd, nbytes, SEEK_CUR);
832  NAMD_write(fd, (char *) &out_integer, sizeof(int32));
833 
834  return(0);
835 }
836 
837 /* A function to update the header. */
838 /* Once all slaves have written the master updates the header. */
839 int update_dcdstep_par_header(int fd)
840 {
841  int32 NSAVC,NSTEP,NFILE;
842  /* don't update header until after write succeeds */
843  OFF_T end = LSEEK(fd,0,SEEK_CUR);
844  LSEEK(fd,NSAVC_POS,SEEK_SET);
845  READ(fd,(void*) &NSAVC,sizeof(int32));
846  LSEEK(fd,NSTEP_POS,SEEK_SET);
847  READ(fd,(void*) &NSTEP,sizeof(int32));
848  LSEEK(fd,NFILE_POS,SEEK_SET);
849  READ(fd,(void*) &NFILE,sizeof(int32));
850  NSTEP += NSAVC;
851  NFILE += 1;
852  LSEEK(fd,NSTEP_POS,SEEK_SET);
853  NAMD_write(fd,(char*) &NSTEP,sizeof(int32));
854  LSEEK(fd,NFILE_POS,SEEK_SET);
855  NAMD_write(fd,(char*) &NFILE,sizeof(int32));
856  LSEEK(fd,end,SEEK_SET);
857 
858  return(0);
859 }
860 
861 /* Each slave writer has a sequential block of atoms to output */
862 /* The stream file position needs to be put at the beginning of each */
863 /* timestep X/Y/Z output */
864 int write_dcdstep_par_slave(int fd, int parL, int parU, int N, float *X, float *Y, float *Z){
865 
866  /* Coordinates for the N elements handled by this writer */
867  int parN = parU-parL+1;
868  size_t nbytes = ((size_t)parN) * 4;
869 
870  /* x's 1st number of Xs */
871  /* skip field for the bytes in X, and the first parL atoms in X*/
872  OFF_T xoffset = sizeof(int)+sizeof(float)*((OFF_T)parL);
873  LSEEK(fd, xoffset, SEEK_CUR);
874  NAMD_write(fd, (char *) X, nbytes);
875 
876  /* skip field for the bytes in X and Y; */
877  /* skip the remaining atoms in X at number of (N-1)-(parU+1)+1
878  * where N-1 is the last atom id, paru+1 is the next atom id. */
879  /* skip the first parL atoms in Y; */
880  OFF_T yoffset = 2*sizeof(int)+sizeof(float)*((OFF_T)(N-parU+parL-1));
881  LSEEK(fd, yoffset, SEEK_CUR);
882  NAMD_write(fd, (char *) Y, nbytes);
883 
884  OFF_T zoffset = yoffset;
885  LSEEK(fd, zoffset, SEEK_CUR);
886  NAMD_write(fd, (char *) Z, nbytes);
887  return(0);
888 }
889 
890 /*****************************************************************************/
891 /* */
892 /* FUNCTION write_dcdheader */
893 /* */
894 /* INPUTS: */
895 /* fd - file descriptor for the dcd file */
896 /* filename - filename for output */
897 /* N - Number of atoms */
898 /* NFILE - Number of sets of coordinates */
899 /* NPRIV - Starting timestep of DCD file - NOT ZERO */
900 /* NSAVC - Timesteps between DCD saves */
901 /* NSTEP - Number of timesteps */
902 /* DELTA - length of a timestep */
903 /* */
904 /* OUTPUTS: */
905 /* none */
906 /* */
907 /* This function prints the "header" information to the DCD file. Since*/
908 /* this is duplicating an unformatted binary output from FORTRAN, its ugly.*/
909 /* So if you're squimish, don't look. */
910 /* Whenever you are updating this function by removing or adding fields, please */
911 /* also update get_dcdheader_size function to reflect the changes! */
912 /* */
913 /*****************************************************************************/
914 
915 int write_dcdheader(int fd, const char *filename, int N, int NFILE, int NPRIV,
916  int NSAVC, int NSTEP, double DELTA, int with_unitcell)
917 {
918  int32 out_integer;
919  float out_float;
920  char title_string[200];
921  int user_id;
922 #ifndef WIN32
923  struct passwd *pwbuf;
924 #endif
925  time_t cur_time;
926  struct tm *tmbuf;
927  char time_str[11];
928 
929  out_integer = 84;
930  NAMD_write(fd, (char *) & out_integer, sizeof(int32));
931  strcpy(title_string, "CORD");
932  NAMD_write(fd, title_string, 4);
933  out_integer = NFILE; /* located at fpos 8 */
934  out_integer = 0; /* ignore the lies */
935  NAMD_write(fd, (char *) & out_integer, sizeof(int32));
936  out_integer = NPRIV; /* located at fpos 12 */
937  NAMD_write(fd, (char *) & out_integer, sizeof(int32));
938  out_integer = NSAVC; /* located at fpos 16 */
939  NAMD_write(fd, (char *) & out_integer, sizeof(int32));
940  out_integer = NSTEP; /* located at fpos 20 */
941  out_integer = NPRIV - NSAVC; /* ignore the lies */
942  NAMD_write(fd, (char *) &out_integer, sizeof(int32));
943  out_integer=0;
944  NAMD_write(fd, (char *) &out_integer, sizeof(int32));
945  NAMD_write(fd, (char *) &out_integer, sizeof(int32));
946  NAMD_write(fd, (char *) &out_integer, sizeof(int32));
947  NAMD_write(fd, (char *) &out_integer, sizeof(int32));
948  NAMD_write(fd, (char *) &out_integer, sizeof(int32));
949  out_float = DELTA;
950  NAMD_write(fd, (char *) &out_float, sizeof(float));
951  out_integer = with_unitcell ? 1 : 0;
952  NAMD_write(fd, (char *) &out_integer, sizeof(int32));
953  out_integer = 0;
954  NAMD_write(fd, (char *) &out_integer, sizeof(int32));
955  NAMD_write(fd, (char *) &out_integer, sizeof(int32));
956  NAMD_write(fd, (char *) &out_integer, sizeof(int32));
957  NAMD_write(fd, (char *) &out_integer, sizeof(int32));
958  NAMD_write(fd, (char *) &out_integer, sizeof(int32));
959  NAMD_write(fd, (char *) &out_integer, sizeof(int32));
960  NAMD_write(fd, (char *) &out_integer, sizeof(int32));
961  NAMD_write(fd, (char *) &out_integer, sizeof(int32));
962  out_integer = 24; // PRETEND TO BE CHARMM24 -JCP
963  NAMD_write(fd, (char *) &out_integer, sizeof(int32));
964  out_integer = 84;
965  NAMD_write(fd, (char *) & out_integer, sizeof(int32));
966 
967  out_integer = 164;
968  NAMD_write(fd, (char *) & out_integer, sizeof(int32));
969  out_integer = 2;
970  NAMD_write(fd, (char *) & out_integer, sizeof(int32));
971 
972  sprintf(title_string, "REMARKS FILENAME=%s CREATED BY NAMD", filename);
973  pad(title_string, 80);
974  NAMD_write(fd, title_string, 80);
975 
976  char username[100];
977 #if defined(WIN32) || defined(NO_GETPWUID)
978  sprintf(username,"Win32");
979 #else
980  user_id= (int) getuid();
981  pwbuf=getpwuid(user_id);
982  if ( pwbuf ) sprintf(username,"%s",pwbuf->pw_name);
983  else sprintf(username,"%d",user_id);
984 #endif
985  cur_time=time(NULL);
986  tmbuf=localtime(&cur_time);
987  strftime(time_str, 10, "%m/%d/%y", tmbuf);
988 
989  sprintf(title_string, "REMARKS DATE: %s CREATED BY USER: %s",
990  time_str, username);
991  pad(title_string, 80);
992  NAMD_write(fd, title_string, 80);
993  out_integer = 164;
994  NAMD_write(fd, (char *) & out_integer, sizeof(int32));
995  out_integer = 4;
996  NAMD_write(fd, (char *) & out_integer, sizeof(int32));
997  out_integer = N;
998  NAMD_write(fd, (char *) & out_integer, sizeof(int32));
999  out_integer = 4;
1000  NAMD_write(fd, (char *) & out_integer, sizeof(int32));
1001 
1002  return(0);
1003 }
1004 
1005 int get_dcdheader_size(){
1006  int headersize = 0;
1007  int totalInt32s = 27; /* 27 writes from out_integer */
1008  int totalChars = 164; /* 3 writes from title_string */
1009  int totalFloats = 1; /* 1 write from out_float */
1010  headersize = sizeof(int32)*totalInt32s+totalChars+sizeof(float)*totalFloats;
1011  return headersize;
1012 }
1013 
1014 /****************************************************************/
1015 /* */
1016 /* FUNCTION close_dcd_read */
1017 /* */
1018 /* INPUTS: */
1019 /* fd - file descriptor to close */
1020 /* num_fixed - the number of fixed atoms */
1021 /* indexes - Array of indexes to be deallocated */
1022 /* */
1023 /* OUTPUTS: */
1024 /* the file pointed to by fd is closed and the memory */
1025 /* pointed to by indexes is freed if any was allocated */
1026 /* */
1027 /* close_dcd_read closes a dcd file that was opened for */
1028 /* reading. It also deallocated memory used for the indexes */
1029 /* used for the free atom list, if there were fixed atoms. */
1030 /* */
1031 /****************************************************************/
1032 
1033 void close_dcd_read(int fd, int num_fixed, int *indexes)
1034 
1035 {
1036 #ifdef WIN32
1037  _close(fd);
1038 #else
1039  close(fd);
1040 #endif
1041 
1042  if (num_fixed)
1043  {
1044  delete [] indexes;
1045  }
1046 }
1047 
1048 /****************************************************************/
1049 /* */
1050 /* FUNCTION close_dcd_write */
1051 /* */
1052 /* INPUTS: */
1053 /* fd - file descriptor to close */
1054 /* */
1055 /* OUTPUTS: */
1056 /* the file pointed to by fd */
1057 /* */
1058 /* close_dcd_write close a dcd file that was opened for */
1059 /* writing */
1060 /* */
1061 /****************************************************************/
1062 
1063 void close_dcd_write(int fd)
1064 
1065 {
1066 #ifdef WIN32
1067  if ( _close(fd) )
1068 #else
1069  if ( fsync(fd) || close(fd) )
1070 #endif
1071  {
1072  NAMD_err("Error closing DCD file");
1073  }
1074 }
1075 
OFF_T
#define OFF_T
Definition: dcdlib.h:42
CHECK_FEOF
#define CHECK_FEOF(X, msg)
Definition: dcdlib.C:97
open_dcd_write_par_slave
int open_dcd_write_par_slave(char *dcdname)
Definition: dcdlib.C:689
open_dcd_read
int open_dcd_read(char *filename)
Definition: dcdlib.C:145
NSAVC_POS
#define NSAVC_POS
Definition: dcdlib.C:630
NAMD_err
void NAMD_err(const char *err_msg)
Definition: common.C:170
update_dcdstep_par_header
int update_dcdstep_par_header(int fd)
Definition: dcdlib.C:839
int32
int32_t int32
Definition: common.h:38
X
#define X
Definition: msm_defn.h:29
NAMD_write
#define NAMD_write
Definition: dcdlib.C:24
CHECK_FREAD
#define CHECK_FREAD(X, msg)
Definition: dcdlib.C:80
O_LARGEFILE
#define O_LARGEFILE
Definition: dcdlib.C:64
open_dcd_write
int open_dcd_write(const char *dcdname)
Definition: dcdlib.C:662
DCD_BADFORMAT
#define DCD_BADFORMAT
Definition: dcdlib.h:51
read_dcdstep
int read_dcdstep(int, int, float *, float *, float *, int, int, int *)
write_dcdstep_par_XYZUnits
int write_dcdstep_par_XYZUnits(int fd, int N)
Definition: dcdlib.C:810
get_dcdheader_size
int get_dcdheader_size()
Definition: dcdlib.C:1005
write_dcdstep_par_cell
int write_dcdstep_par_cell(int fd, double *cell)
Definition: dcdlib.C:785
NAMD_seek
OFF_T NAMD_seek(int file, OFF_T offset, int whence)
Definition: dcdlib.C:49
NAMD_bug
void NAMD_bug(const char *err_msg)
Definition: common.C:195
Z
#define Z
Definition: msm_defn.h:33
close_dcd_read
void close_dcd_read(int fd, int num_fixed, int *indexes)
Definition: dcdlib.C:1033
READ
#define READ
Definition: dcdlib.C:46
write_dcdstep
int write_dcdstep(int fd, int N, float *X, float *Y, float *Z, double *cell)
Definition: dcdlib.C:736
NAMD_die
void NAMD_die(const char *err_msg)
Definition: common.C:147
NFILE_POS
#define NFILE_POS
Definition: dcdlib.C:628
write_dcdheader
int write_dcdheader(int fd, const char *filename, int N, int NFILE, int NPRIV, int NSAVC, int NSTEP, double DELTA, int with_unitcell)
Definition: dcdlib.C:915
NAMD_backup_file
void NAMD_backup_file(const char *filename, const char *extension)
Definition: common.C:235
NSTEP_POS
#define NSTEP_POS
Definition: dcdlib.C:631
DCD_DNE
#define DCD_DNE
Definition: dcdlib.h:47
write_dcdstep_par_slave
int write_dcdstep_par_slave(int fd, int parL, int parU, int N, float *X, float *Y, float *Z)
Definition: dcdlib.C:864
Y
#define Y
Definition: msm_defn.h:31
close_dcd_write
void close_dcd_write(int fd)
Definition: dcdlib.C:1063
pad
void pad(char *s, int len)
Definition: dcdlib.C:105
DCD_BADMALLOC
#define DCD_BADMALLOC
Definition: dcdlib.h:53
read_dcdheader
int read_dcdheader(int, int *, int *, int *, int *, double *, int *, int **)
DCD_OPENFAILED
#define DCD_OPENFAILED
Definition: dcdlib.h:48
LSEEK
#define LSEEK
Definition: dcdlib.C:61
Generated by   doxygen 1.8.14