Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : University of Illinois Open Source License
3 : Copyright 2003 Theoretical and Computational Biophysics Group,
4 : All rights reserved.
5 :
6 : Developed by: Theoretical and Computational Biophysics Group
7 : University of Illinois at Urbana-Champaign
8 : http://www.ks.uiuc.edu/
9 :
10 : Permission is hereby granted, free of charge, to any person obtaining a copy of
11 : this software and associated documentation files (the Software), to deal with
12 : the Software without restriction, including without limitation the rights to
13 : use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
14 : of the Software, and to permit persons to whom the Software is furnished to
15 : do so, subject to the following conditions:
16 :
17 : Redistributions of source code must retain the above copyright notice,
18 : this list of conditions and the following disclaimers.
19 :
20 : Redistributions in binary form must reproduce the above copyright notice,
21 : this list of conditions and the following disclaimers in the documentation
22 : and/or other materials provided with the distribution.
23 :
24 : Neither the names of Theoretical and Computational Biophysics Group,
25 : University of Illinois at Urbana-Champaign, nor the names of its contributors
26 : may be used to endorse or promote products derived from this Software without
27 : specific prior written permission.
28 :
29 : THE SOFTWARE IS PROVIDED AS IS, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
30 : IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
31 : FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
32 : THE CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
33 : OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
34 : ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
35 : OTHER DEALINGS WITH THE SOFTWARE.
36 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
37 : #if defined(__PLUMED_HAS_MOLFILE_PLUGINS) && ! defined(__PLUMED_HAS_EXTERNAL_MOLFILE_PLUGINS)
38 : /***************************************************************************
39 : *cr
40 : *cr (C) Copyright 1995-2016 The Board of Trustees of the
41 : *cr University of Illinois
42 : *cr All Rights Reserved
43 : *cr
44 : ***************************************************************************/
45 :
46 : /***************************************************************************
47 : * RCS INFORMATION:
48 : *
49 : * $RCSfile: dcdplugin.c,v $
50 : * $Author: johns $ $Locker: $ $State: Exp $
51 : * $Revision: 1.88 $ $Date: 2020/12/17 17:14:07 $
52 : *
53 : ***************************************************************************
54 : * DESCRIPTION:
55 : * Code for reading and writing CHARMM, NAMD, and X-PLOR format
56 : * molecular dynamic trajectory files.
57 : *
58 : * TODO:
59 : * Integrate improvements from the NAMD source tree
60 : * - NAMD's writer code has better type-correctness for the sizes
61 : * of "int". NAMD uses "int32" explicitly, which is required on
62 : * machines like the T3E. VMD's version of the code doesn't do that
63 : * presently.
64 : *
65 : * Try various alternative I/O API options:
66 : * - use mmap(), with read-once flags
67 : * - use O_DIRECT open mode on new revs of Linux kernel
68 : * - use directio() call on a file descriptor to enable on Solaris
69 : * - use aio_open()/read()/write()
70 : * - use readv/writev() etc.
71 : *
72 : * Standalone test binary compilation flags:
73 : * cc -fast -xarch=v9a -I../../include -DTEST_DCDPLUGIN dcdplugin.c \
74 : * -o ~/bin/readdcd -lm
75 : *
76 : * Profiling flags:
77 : * cc -xpg -fast -xarch=v9a -g -I../../include -DTEST_DCDPLUGIN dcdplugin.c \
78 : * -o ~/bin/readdcd -lm
79 : *
80 : ***************************************************************************/
81 :
82 : #include "largefiles.h" /* platform dependent 64-bit file I/O defines */
83 : #include "fastio.h" /* must come before others, for O_DIRECT... */
84 :
85 : #include <stdio.h>
86 : #include <sys/stat.h>
87 : #include <sys/types.h>
88 : #include <stdlib.h>
89 : #include <stddef.h>
90 : #include <string.h>
91 : #include <math.h>
92 : #include <time.h>
93 : #include "endianswap.h"
94 : #include "molfile_plugin.h"
95 :
96 : namespace PLMD{
97 : namespace molfile{
98 :
99 : #ifndef M_PI_2
100 : #define M_PI_2 1.57079632679489661922
101 : #endif
102 :
103 : #define RECSCALE32BIT 1
104 : #define RECSCALE64BIT 2
105 : #define RECSCALEMAX 2
106 :
107 : typedef struct {
108 : fio_fd fd;
109 : int natoms;
110 : int nsets;
111 : int setsread;
112 : int istart;
113 : int nsavc;
114 : double delta;
115 : int nfixed;
116 : float *x, *y, *z;
117 : int *freeind;
118 : float *fixedcoords;
119 : int reverse;
120 : int charmm;
121 : int first;
122 : int with_unitcell;
123 : } dcdhandle;
124 :
125 : /* Define error codes that may be returned by the DCD routines */
126 : #define DCD_SUCCESS 0 /* No problems */
127 : #define DCD_EOF -1 /* Normal EOF */
128 : #define DCD_DNE -2 /* DCD file does not exist */
129 : #define DCD_OPENFAILED -3 /* Open of DCD file failed */
130 : #define DCD_BADREAD -4 /* read call on DCD file failed */
131 : #define DCD_BADEOF -5 /* premature EOF found in DCD file */
132 : #define DCD_BADFORMAT -6 /* format of DCD file is wrong */
133 : #define DCD_FILEEXISTS -7 /* output file already exists */
134 : #define DCD_BADMALLOC -8 /* malloc failed */
135 : #define DCD_BADWRITE -9 /* write call on DCD file failed */
136 :
137 : /* Define feature flags for this DCD file */
138 : #define DCD_IS_XPLOR 0x00
139 : #define DCD_IS_CHARMM 0x01
140 : #define DCD_HAS_4DIMS 0x02
141 : #define DCD_HAS_EXTRA_BLOCK 0x04
142 : #define DCD_HAS_64BIT_REC 0x08
143 :
144 : /* defines used by write_dcdstep */
145 : #define NFILE_POS 8L
146 : #define NSTEP_POS 20L
147 :
148 : /* READ Macro to make porting easier */
149 : #define READ(fd, buf, size) fio_fread(((void *) buf), (size), 1, (fd))
150 :
151 : /* WRITE Macro to make porting easier */
152 : #define WRITE(fd, buf, size) fio_fwrite(((void *) buf), (size), 1, (fd))
153 :
154 : /* XXX This is broken - fread never returns -1 */
155 : #define CHECK_FREAD(X, msg) if (X==-1) { return(DCD_BADREAD); }
156 : #define CHECK_FEOF(X, msg) if (X==0) { return(DCD_BADEOF); }
157 :
158 :
159 : /* print DCD error in a human readable way */
160 0 : static void print_dcderror(const char *func, int errcode) {
161 : const char *errstr;
162 :
163 0 : switch (errcode) {
164 : case DCD_EOF: errstr = "end of file"; break;
165 0 : case DCD_DNE: errstr = "file not found"; break;
166 0 : case DCD_OPENFAILED: errstr = "file open failed"; break;
167 0 : case DCD_BADREAD: errstr = "error during read"; break;
168 0 : case DCD_BADEOF: errstr = "premature end of file"; break;
169 0 : case DCD_BADFORMAT: errstr = "corruption or unrecognized file structure"; break;
170 0 : case DCD_FILEEXISTS: errstr = "output file already exists"; break;
171 0 : case DCD_BADMALLOC: errstr = "memory allocation failed"; break;
172 0 : case DCD_BADWRITE: errstr = "error during write"; break;
173 0 : case DCD_SUCCESS:
174 : default:
175 : errstr = "no error";
176 0 : break;
177 : }
178 : printf("dcdplugin) %s: %s\n", func, errstr);
179 0 : }
180 :
181 :
182 : /*
183 : * Read the header information from a dcd file.
184 : * Input: fd - a file struct opened for binary reading.
185 : * Output: 0 on success, negative error code on failure.
186 : * Side effects: *natoms set to number of atoms per frame
187 : * *nsets set to number of frames in dcd file
188 : * *istart set to starting timestep of dcd file
189 : * *nsavc set to timesteps between dcd saves
190 : * *delta set to value of trajectory timestep
191 : * *nfixed set to number of fixed atoms
192 : * *freeind may be set to heap-allocated space
193 : * *reverse set to one if reverse-endian, zero if not.
194 : * *charmm set to internal code for handling charmm data.
195 : */
196 1 : static int read_dcdheader(fio_fd fd, int *N, int *NSET, int *ISTART,
197 : int *NSAVC, double *DELTA, int *NAMNF,
198 : int **FREEINDEXES, float **fixedcoords, int *reverseEndian,
199 : int *charmm)
200 : {
201 : unsigned int input_integer[2]; /* buffer space */
202 : int i, ret_val, rec_scale;
203 : char hdrbuf[84]; /* char buffer used to store header */
204 : int NTITLE;
205 : int dcdcordmagic;
206 : int hugefile = 0;
207 : char *corp = (char *) &dcdcordmagic;
208 :
209 : /* coordinate dcd file magic string 'CORD' */
210 : corp[0] = 'C';
211 : corp[1] = 'O';
212 : corp[2] = 'R';
213 : corp[3] = 'D';
214 :
215 : /* First thing in the file should be an 84.
216 : * some 64-bit compiles have a 64-bit record length indicator,
217 : * so we have to read two ints and check in a more complicated
218 : * way. :-( */
219 1 : ret_val = READ(fd, input_integer, 2*sizeof(unsigned int));
220 1 : CHECK_FREAD(ret_val, "reading first int from dcd file");
221 1 : CHECK_FEOF(ret_val, "reading first int from dcd file");
222 :
223 : /* Check magic number in file header and determine byte order*/
224 1 : if ((input_integer[0]+input_integer[1]) == 84) {
225 0 : *reverseEndian=0;
226 : rec_scale=RECSCALE64BIT;
227 : printf("dcdplugin) detected CHARMM -i8 64-bit DCD file of native endianness\n");
228 1 : } else if (input_integer[0] == 84 && input_integer[1] == dcdcordmagic) {
229 1 : *reverseEndian=0;
230 : rec_scale=RECSCALE32BIT;
231 : printf("dcdplugin) detected standard 32-bit DCD file of native endianness\n");
232 : } else {
233 : /* now try reverse endian */
234 0 : swap4_aligned(input_integer, 2); /* will have to unswap magic if 32-bit */
235 0 : if ((input_integer[0]+input_integer[1]) == 84) {
236 0 : *reverseEndian=1;
237 : rec_scale=RECSCALE64BIT;
238 : printf("dcdplugin) detected CHARMM -i8 64-bit DCD file of opposite endianness\n");
239 : } else {
240 0 : swap4_aligned(&input_integer[1], 1); /* unswap magic (see above) */
241 0 : if (input_integer[0] == 84 && input_integer[1] == dcdcordmagic) {
242 0 : *reverseEndian=1;
243 : rec_scale=RECSCALE32BIT;
244 : printf("dcdplugin) detected standard 32-bit DCD file of opposite endianness\n");
245 : } else {
246 : /* not simply reversed endianism or -i8, something rather more evil */
247 : printf("dcdplugin) unrecognized DCD header:\n");
248 0 : printf("dcdplugin) [0]: %10d [1]: %10d\n", input_integer[0], input_integer[1]);
249 0 : printf("dcdplugin) [0]: 0x%08x [1]: 0x%08x\n", input_integer[0], input_integer[1]);
250 0 : return DCD_BADFORMAT;
251 :
252 : }
253 : }
254 : }
255 :
256 : /* check for magic string, in case of long record markers */
257 : if (rec_scale == RECSCALE64BIT) {
258 0 : ret_val = READ(fd, input_integer, sizeof(unsigned int));
259 0 : if (input_integer[0] != dcdcordmagic) {
260 : printf("dcdplugin) failed to find CORD magic in CHARMM -i8 64-bit DCD file\n");
261 0 : return DCD_BADFORMAT;
262 : }
263 : }
264 :
265 : /* Buffer the entire header for random access */
266 1 : ret_val = READ(fd, hdrbuf, 80);
267 1 : CHECK_FREAD(ret_val, "buffering header");
268 1 : CHECK_FEOF(ret_val, "buffering header");
269 :
270 : /* CHARMm-genereate DCD files set the last integer in the */
271 : /* header, which is unused by X-PLOR, to its version number. */
272 : /* Checking if this is nonzero tells us this is a CHARMm file */
273 : /* and to look for other CHARMm flags. */
274 1 : if (*((int *) (hdrbuf + 76)) != 0) {
275 1 : (*charmm) = DCD_IS_CHARMM;
276 1 : if (*((int *) (hdrbuf + 40)) != 0)
277 1 : (*charmm) |= DCD_HAS_EXTRA_BLOCK;
278 :
279 1 : if (*((int *) (hdrbuf + 44)) == 1)
280 0 : (*charmm) |= DCD_HAS_4DIMS;
281 :
282 1 : if (rec_scale == RECSCALE64BIT)
283 0 : (*charmm) |= DCD_HAS_64BIT_REC;
284 :
285 : } else {
286 0 : (*charmm) = DCD_IS_XPLOR; /* must be an X-PLOR format DCD file */
287 : }
288 :
289 1 : if (*charmm & DCD_IS_CHARMM) {
290 : /* CHARMM and NAMD versions 2.1b1 and later */
291 : printf("dcdplugin) CHARMM format DCD file (also NAMD 2.1 and later)\n");
292 : } else {
293 : /* CHARMM and NAMD versions prior to 2.1b1 */
294 : printf("dcdplugin) X-PLOR format DCD file (also NAMD 2.0 and earlier)\n");
295 : }
296 :
297 : /* Store the number of sets of coordinates (NSET) */
298 1 : (*NSET) = *((int *) (hdrbuf));
299 1 : if (*reverseEndian) swap4_unaligned(NSET, 1);
300 :
301 : /* Store ISTART, the starting timestep */
302 1 : (*ISTART) = *((int *) (hdrbuf + 4));
303 1 : if (*reverseEndian) swap4_unaligned(ISTART, 1);
304 :
305 : /* Store NSAVC, the number of timesteps between dcd saves */
306 1 : (*NSAVC) = *((int *) (hdrbuf + 8));
307 1 : if (*reverseEndian) swap4_unaligned(NSAVC, 1);
308 :
309 : /* Store NAMNF, the number of fixed atoms */
310 1 : (*NAMNF) = *((int *) (hdrbuf + 32));
311 1 : if (*reverseEndian) swap4_unaligned(NAMNF, 1);
312 :
313 : /* Read in the timestep, DELTA */
314 : /* Note: DELTA is stored as a double with X-PLOR but as a float with CHARMm */
315 1 : if ((*charmm) & DCD_IS_CHARMM) {
316 : float ftmp;
317 1 : ftmp = *((float *)(hdrbuf+36)); /* is this safe on Alpha? */
318 1 : if (*reverseEndian)
319 0 : swap4_aligned(&ftmp, 1);
320 :
321 1 : *DELTA = (double)ftmp;
322 : } else {
323 0 : (*DELTA) = *((double *)(hdrbuf + 36));
324 0 : if (*reverseEndian) swap8_unaligned(DELTA, 1);
325 : }
326 :
327 : /* Get the end size of the first block */
328 1 : ret_val = READ(fd, input_integer, rec_scale*sizeof(int));
329 1 : CHECK_FREAD(ret_val, "reading second 84 from dcd file");
330 1 : CHECK_FEOF(ret_val, "reading second 84 from dcd file");
331 1 : if (*reverseEndian) swap4_aligned(input_integer, rec_scale);
332 :
333 1 : if (rec_scale == RECSCALE64BIT) {
334 0 : if ((input_integer[0]+input_integer[1]) != 84) {
335 : return DCD_BADFORMAT;
336 : }
337 : } else {
338 1 : if (input_integer[0] != 84) {
339 : return DCD_BADFORMAT;
340 : }
341 : }
342 :
343 : /* Read in the size of the next block */
344 1 : input_integer[1] = 0;
345 1 : ret_val = READ(fd, input_integer, rec_scale*sizeof(int));
346 1 : CHECK_FREAD(ret_val, "reading size of title block");
347 1 : CHECK_FEOF(ret_val, "reading size of title block");
348 1 : if (*reverseEndian) swap4_aligned(input_integer, rec_scale);
349 :
350 1 : if ((((input_integer[0]+input_integer[1])-4) % 80) == 0) {
351 : /* Read NTITLE, the number of 80 character title strings there are */
352 1 : ret_val = READ(fd, &NTITLE, sizeof(int));
353 1 : CHECK_FREAD(ret_val, "reading NTITLE");
354 1 : CHECK_FEOF(ret_val, "reading NTITLE");
355 1 : if (*reverseEndian) swap4_aligned(&NTITLE, 1);
356 :
357 1 : if (NTITLE < 0) {
358 : printf("dcdplugin) WARNING: Bogus NTITLE value: %d (hex: %08x)\n",
359 : NTITLE, NTITLE);
360 0 : return DCD_BADFORMAT;
361 : }
362 :
363 1 : if (NTITLE > 1000) {
364 : printf("dcdplugin) WARNING: Bogus NTITLE value: %d (hex: %08x)\n",
365 : NTITLE, NTITLE);
366 0 : if (NTITLE == 1095062083) {
367 : printf("dcdplugin) WARNING: Broken Vega ZZ 2.4.0 DCD file detected\n");
368 : printf("dcdplugin) Assuming 2 title lines, good luck...\n");
369 0 : NTITLE = 2;
370 : } else {
371 : printf("dcdplugin) Assuming zero title lines, good luck...\n");
372 0 : NTITLE = 0;
373 : }
374 : }
375 :
376 3 : for (i=0; i<NTITLE; i++) {
377 2 : fio_fseek(fd, 80, FIO_SEEK_CUR);
378 : CHECK_FEOF(ret_val, "reading TITLE");
379 : }
380 :
381 : /* Get the ending size for this block */
382 1 : ret_val = READ(fd, input_integer, rec_scale*sizeof(int));
383 1 : CHECK_FREAD(ret_val, "reading size of title block");
384 1 : CHECK_FEOF(ret_val, "reading size of title block");
385 : } else {
386 : return DCD_BADFORMAT;
387 : }
388 :
389 : /* Read in an integer '4' */
390 1 : input_integer[1] = 0;
391 1 : ret_val = READ(fd, input_integer, rec_scale*sizeof(int));
392 :
393 1 : CHECK_FREAD(ret_val, "reading a '4'");
394 1 : CHECK_FEOF(ret_val, "reading a '4'");
395 1 : if (*reverseEndian) swap4_aligned(input_integer, rec_scale);
396 :
397 1 : if ((input_integer[0]+input_integer[1]) != 4) {
398 : return DCD_BADFORMAT;
399 : }
400 :
401 : /* Read in the number of atoms */
402 1 : ret_val = READ(fd, N, sizeof(int));
403 1 : CHECK_FREAD(ret_val, "reading number of atoms");
404 1 : CHECK_FEOF(ret_val, "reading number of atoms");
405 1 : if (*reverseEndian) swap4_aligned(N, 1);
406 :
407 1 : if (*N > (1L<<30)) {
408 : hugefile=1;
409 : printf("dcdplugin) ***\n");
410 : printf("dcdplugin) *** Trajectory contains over 2^30 atoms.\n");
411 : printf("dcdplugin) *** Huge file integer wraparound handling enabled.\n");
412 : printf("dcdplugin) ***\n");
413 : }
414 :
415 : /* Read in an integer '4' */
416 1 : input_integer[1] = 0;
417 1 : ret_val = READ(fd, input_integer, rec_scale*sizeof(int));
418 1 : CHECK_FREAD(ret_val, "reading a '4'");
419 1 : CHECK_FEOF(ret_val, "reading a '4'");
420 1 : if (*reverseEndian) swap4_aligned(input_integer, rec_scale);
421 :
422 1 : if ((input_integer[0]+input_integer[1]) != 4) {
423 : return DCD_BADFORMAT;
424 : }
425 :
426 1 : *FREEINDEXES = NULL;
427 1 : *fixedcoords = NULL;
428 1 : if (*NAMNF != 0) {
429 0 : (*FREEINDEXES) = (int *) calloc((((ptrdiff_t)(*N))-(*NAMNF)), sizeof(int));
430 0 : if (*FREEINDEXES == NULL)
431 : return DCD_BADMALLOC;
432 :
433 0 : *fixedcoords = (float *) calloc(((ptrdiff_t)(*N))*4L - (*NAMNF), sizeof(float));
434 0 : if (*fixedcoords == NULL)
435 : return DCD_BADMALLOC;
436 :
437 : /* Read in index array size */
438 0 : input_integer[1]=0;
439 0 : ret_val = READ(fd, input_integer, rec_scale*sizeof(int));
440 0 : CHECK_FREAD(ret_val, "reading size of index array");
441 0 : CHECK_FEOF(ret_val, "reading size of index array");
442 0 : if (*reverseEndian) swap4_aligned(input_integer, rec_scale);
443 :
444 : /* when we have more then 2^30 atoms, tests like this one */
445 : /* are no longer meaningful and have to be bypassed... */
446 0 : if (!hugefile && ((input_integer[0]+input_integer[1]) != ((*N)-(*NAMNF))*4L)) {
447 : return DCD_BADFORMAT;
448 : }
449 :
450 0 : ret_val = READ(fd, (*FREEINDEXES), ((ptrdiff_t) ((*N)-(*NAMNF)))*sizeof(int));
451 0 : CHECK_FREAD(ret_val, "reading size of index array");
452 0 : CHECK_FEOF(ret_val, "reading size of index array");
453 :
454 0 : if (*reverseEndian)
455 0 : swap4_aligned((*FREEINDEXES), ((*N)-(*NAMNF)));
456 :
457 0 : input_integer[1]=0;
458 0 : ret_val = READ(fd, input_integer, rec_scale*sizeof(int));
459 0 : CHECK_FREAD(ret_val, "reading size of index array");
460 0 : CHECK_FEOF(ret_val, "reading size of index array");
461 0 : if (*reverseEndian) swap4_aligned(input_integer, rec_scale);
462 :
463 : /* when we have more then 2^30 atoms, tests like this one */
464 : /* are no longer meaningful and have to be bypassed... */
465 0 : if (!hugefile && ((input_integer[0]+input_integer[1]) != ((*N)-(*NAMNF))*4L)) {
466 0 : return DCD_BADFORMAT;
467 : }
468 : }
469 :
470 : return DCD_SUCCESS;
471 : }
472 :
473 :
474 21 : static int read_charmm_extrablock(fio_fd fd, int charmm, int reverseEndian,
475 : float *unitcell) {
476 : int i, input_integer[2], rec_scale;
477 :
478 21 : if (charmm & DCD_HAS_64BIT_REC) {
479 : rec_scale = RECSCALE64BIT;
480 : } else {
481 : rec_scale = RECSCALE32BIT;
482 : }
483 :
484 21 : if ((charmm & DCD_IS_CHARMM) && (charmm & DCD_HAS_EXTRA_BLOCK)) {
485 : /* Leading integer must be 48 */
486 21 : input_integer[1] = 0;
487 21 : if (fio_fread(input_integer, sizeof(int), rec_scale, fd) != rec_scale)
488 : return DCD_BADREAD;
489 21 : if (reverseEndian) swap4_aligned(input_integer, rec_scale);
490 21 : if ((input_integer[0]+input_integer[1]) == 48) {
491 : double tmp[6];
492 21 : if (fio_fread(tmp, 48, 1, fd) != 1) return DCD_BADREAD;
493 21 : if (reverseEndian)
494 0 : swap8_aligned(tmp, 6);
495 147 : for (i=0; i<6; i++) unitcell[i] = (float)tmp[i];
496 : } else {
497 : /* unrecognized block, just skip it */
498 0 : if (fio_fseek(fd, (input_integer[0]+input_integer[1]), FIO_SEEK_CUR)) return DCD_BADREAD;
499 : }
500 21 : if (fio_fread(input_integer, sizeof(int), rec_scale, fd) != rec_scale) return DCD_BADREAD;
501 : }
502 :
503 : return DCD_SUCCESS;
504 : }
505 :
506 :
507 0 : static int read_fixed_atoms(fio_fd fd, int N, int num_free, const int *indexes,
508 : int reverseEndian, const float *fixedcoords,
509 : float *freeatoms, float *pos, int charmm) {
510 : int i, input_integer[2], rec_scale;
511 :
512 0 : if(charmm & DCD_HAS_64BIT_REC) {
513 : rec_scale=RECSCALE64BIT;
514 : } else {
515 : rec_scale=RECSCALE32BIT;
516 : }
517 :
518 : /* Read leading integer */
519 0 : input_integer[1]=0;
520 0 : if (fio_fread(input_integer, sizeof(int), rec_scale, fd) != rec_scale) return DCD_BADREAD;
521 0 : if (reverseEndian) swap4_aligned(input_integer, rec_scale);
522 0 : if ((input_integer[0]+input_integer[1]) != 4L*num_free) return DCD_BADFORMAT;
523 :
524 : /* Read free atom coordinates */
525 0 : if (fio_fread(freeatoms, 4L*num_free, 1, fd) != 1) return DCD_BADREAD;
526 0 : if (reverseEndian)
527 0 : swap4_aligned(freeatoms, num_free);
528 :
529 : /* Copy fixed and free atom coordinates into position buffer */
530 0 : memcpy(pos, fixedcoords, 4L*N);
531 0 : for (i=0; i<num_free; i++)
532 0 : pos[indexes[i]-1] = freeatoms[i];
533 :
534 : /* Read trailing integer */
535 0 : input_integer[1]=0;
536 0 : if (fio_fread(input_integer, sizeof(int), rec_scale, fd) != rec_scale) return DCD_BADREAD;
537 0 : if (reverseEndian) swap4_aligned(input_integer, rec_scale);
538 0 : if ((input_integer[0]+input_integer[1]) != 4L*num_free) return DCD_BADFORMAT;
539 :
540 : return DCD_SUCCESS;
541 : }
542 :
543 :
544 21 : static int read_charmm_4dim(fio_fd fd, int charmm, int reverseEndian) {
545 : int input_integer[2], rec_scale;
546 :
547 21 : if (charmm & DCD_HAS_64BIT_REC) {
548 : rec_scale=RECSCALE64BIT;
549 : } else {
550 : rec_scale=RECSCALE32BIT;
551 : }
552 :
553 : /* If this is a CHARMm file and contains a 4th dimension block, */
554 : /* we must skip past it to avoid problems */
555 21 : if ((charmm & DCD_IS_CHARMM) && (charmm & DCD_HAS_4DIMS)) {
556 0 : input_integer[1]=0;
557 0 : if (fio_fread(input_integer, sizeof(int), rec_scale, fd) != rec_scale) return DCD_BADREAD;
558 0 : if (reverseEndian) swap4_aligned(input_integer, rec_scale);
559 0 : if (fio_fseek(fd, (input_integer[0]+input_integer[1]), FIO_SEEK_CUR)) return DCD_BADREAD;
560 0 : if (fio_fread(input_integer, sizeof(int), rec_scale, fd) != rec_scale) return DCD_BADREAD;
561 : }
562 :
563 : return DCD_SUCCESS;
564 : }
565 :
566 :
567 : /*
568 : * Read a dcd timestep from a dcd file
569 : * Input: fd - a file struct opened for binary reading, from which the
570 : * header information has already been read.
571 : * natoms, nfixed, first, *freeind, reverse, charmm - the corresponding
572 : * items as set by read_dcdheader
573 : * first - true if this is the first frame we are reading.
574 : * x, y, z: space for natoms each of floats.
575 : * unitcell - space for six floats to hold the unit cell data.
576 : * Not set if no unit cell data is present.
577 : * Output: 0 on success, negative error code on failure.
578 : * Side effects: x, y, z contain the coordinates for the timestep read.
579 : * unitcell holds unit cell data if present.
580 : */
581 21 : static int read_dcdstep(fio_fd fd, int N, float *X, float *Y, float *Z,
582 : float *unitcell, int num_fixed,
583 : int first, int *indexes, float *fixedcoords,
584 : int reverseEndian, int charmm) {
585 : int ret_val; /* Return value from read */
586 : ptrdiff_t rec_scale;
587 21 : int hugefile = (N > (1L<<30)) ? 1 : 0;
588 : int check_reclen = 1; /* Enable Fortran record length value safety checks */
589 :
590 : /* Fortran record length checks disabled for huge files or user request */
591 21 : if (hugefile || (getenv("VMDDCDNOCHECKRECLEN") != NULL))
592 : check_reclen = 0;
593 :
594 21 : if (charmm & DCD_HAS_64BIT_REC) {
595 : rec_scale=RECSCALE64BIT;
596 : } else {
597 : rec_scale=RECSCALE32BIT;
598 : }
599 :
600 21 : if ((num_fixed==0) || first) {
601 : /* temp storage for reading formatting info */
602 : /* note: has to be max size we'll ever use */
603 : int tmpbuf[6L*RECSCALEMAX];
604 :
605 : fio_iovec iov[7]; /* I/O vector for fio_readv() call */
606 : int i;
607 :
608 : /* if there are no fixed atoms or this is the first timestep read */
609 : /* then we read all coordinates normally. */
610 :
611 : /* read the charmm periodic cell information */
612 : /* XXX this too should be read together with the other items in a */
613 : /* single fio_readv() call in order to prevent lots of extra */
614 : /* kernel/user context switches. */
615 21 : ret_val = read_charmm_extrablock(fd, charmm, reverseEndian, unitcell);
616 21 : if (ret_val) return ret_val;
617 :
618 : /* setup the I/O vector for the call to fio_readv() */
619 : iov[0].iov_base = (fio_caddr_t) &tmpbuf[0]; /* read format integer */
620 21 : iov[0].iov_len = rec_scale*sizeof(int);
621 :
622 : iov[1].iov_base = (fio_caddr_t) X; /* read X coordinates */
623 21 : iov[1].iov_len = sizeof(float)*N;
624 :
625 21 : iov[2].iov_base = (fio_caddr_t) &tmpbuf[1*rec_scale]; /* read 2 format integers */
626 21 : iov[2].iov_len = rec_scale*sizeof(int) * 2L;
627 :
628 : iov[3].iov_base = (fio_caddr_t) Y; /* read Y coordinates */
629 : iov[3].iov_len = sizeof(float)*N;
630 :
631 21 : iov[4].iov_base = (fio_caddr_t) &tmpbuf[3L*rec_scale]; /* read 2 format integers */
632 : iov[4].iov_len = rec_scale*sizeof(int) * 2L;
633 :
634 : iov[5].iov_base = (fio_caddr_t) Z; /* read Y coordinates */
635 : iov[5].iov_len = sizeof(float)*N;
636 :
637 21 : iov[6].iov_base = (fio_caddr_t) &tmpbuf[5L*rec_scale]; /* read format integer */
638 : iov[6].iov_len = rec_scale*sizeof(int);
639 :
640 : #if 1
641 : /* Use fall-back code instead of readv(): */
642 : /* Some platforms implement readv() as user level code in libc, */
643 : /* and due to POSIX atomicity requirements for readv()/writev(), */
644 : /* they may copy data to internal temp buffers, which can kill */
645 : /* performance, and in cases when doing single I/O ops on large, */
646 : /* buffers, e.g. > 2GB, can fail with shorts reads or writes... */
647 : /* On such platforms it is best to avoid using readv()/writev()... */
648 : {
649 : int readcnt = 0;
650 21 : readcnt = fio_fread(iov[0].iov_base, iov[0].iov_len, 1, fd);
651 21 : readcnt += fio_fread(iov[1].iov_base, iov[1].iov_len, 1, fd);
652 21 : readcnt += fio_fread(iov[2].iov_base, iov[2].iov_len, 1, fd);
653 21 : readcnt += fio_fread(iov[3].iov_base, iov[3].iov_len, 1, fd);
654 21 : readcnt += fio_fread(iov[4].iov_base, iov[4].iov_len, 1, fd);
655 21 : readcnt += fio_fread(iov[5].iov_base, iov[5].iov_len, 1, fd);
656 21 : readcnt += fio_fread(iov[6].iov_base, iov[6].iov_len, 1, fd);
657 :
658 : /* if all records read correctly, then the reads are okay */
659 21 : if (readcnt != 7)
660 : return DCD_BADREAD;
661 : }
662 : #else
663 : /* check number of bytes actually read */
664 : if (fio_readv(fd, &iov[0], 7) != ((fio_size_t) (rec_scale*6L*sizeof(int) + 3L*N*sizeof(float))))
665 : return DCD_BADREAD;
666 : #endif
667 :
668 : /* convert endianism if necessary */
669 21 : if (reverseEndian) {
670 0 : swap4_aligned(&tmpbuf[0], rec_scale*6L);
671 0 : swap4_aligned(X, N);
672 0 : swap4_aligned(Y, N);
673 0 : swap4_aligned(Z, N);
674 : }
675 :
676 : /* when we have more then 2^30 atoms, tests like this one */
677 : /* are no longer meaningful and have to be bypassed... */
678 21 : if (check_reclen) {
679 : /* double-check the fortran format size values for safety */
680 21 : if (rec_scale == 1) {
681 147 : for (i=0; i<6; i++) {
682 126 : if (tmpbuf[i] != sizeof(float)*N) return DCD_BADFORMAT;
683 : }
684 : } else {
685 0 : for (i=0; i<6; i++) {
686 0 : if ((tmpbuf[2L*i]+tmpbuf[2L*i+1L]) != sizeof(float)*N) return DCD_BADFORMAT;
687 : }
688 : }
689 : }
690 :
691 : /* copy fixed atom coordinates into fixedcoords array if this was the */
692 : /* first timestep, to be used from now on. We just copy all atoms. */
693 21 : if (num_fixed && first) {
694 : memcpy(fixedcoords, X, N*sizeof(float));
695 0 : memcpy(fixedcoords+N, Y, N*sizeof(float));
696 0 : memcpy(fixedcoords+2L*N, Z, N*sizeof(float));
697 : }
698 :
699 : /* read in the optional charmm 4th array */
700 : /* XXX this too should be read together with the other items in a */
701 : /* single fio_readv() call in order to prevent lots of extra */
702 : /* kernel/user context switches. */
703 21 : ret_val = read_charmm_4dim(fd, charmm, reverseEndian);
704 21 : if (ret_val) return ret_val;
705 : } else {
706 : /* if there are fixed atoms, and this isn't the first frame, then we */
707 : /* only read in the non-fixed atoms for all subsequent timesteps. */
708 0 : ret_val = read_charmm_extrablock(fd, charmm, reverseEndian, unitcell);
709 0 : if (ret_val) return ret_val;
710 0 : ret_val = read_fixed_atoms(fd, N, N-num_fixed, indexes, reverseEndian,
711 0 : fixedcoords, fixedcoords+3L*N, X, charmm);
712 0 : if (ret_val) return ret_val;
713 0 : ret_val = read_fixed_atoms(fd, N, N-num_fixed, indexes, reverseEndian,
714 0 : fixedcoords+N, fixedcoords+3L*N, Y, charmm);
715 0 : if (ret_val) return ret_val;
716 0 : ret_val = read_fixed_atoms(fd, N, N-num_fixed, indexes, reverseEndian,
717 0 : fixedcoords+2*N, fixedcoords+3L*N, Z, charmm);
718 0 : if (ret_val) return ret_val;
719 0 : ret_val = read_charmm_4dim(fd, charmm, reverseEndian);
720 0 : if (ret_val) return ret_val;
721 : }
722 :
723 : return DCD_SUCCESS;
724 : }
725 :
726 :
727 : /*
728 : * Skip past a timestep. If there are fixed atoms, this cannot be used with
729 : * the first timestep.
730 : * Input: fd - a file struct from which the header has already been read
731 : * natoms - number of atoms per timestep
732 : * nfixed - number of fixed atoms
733 : * charmm - charmm flags as returned by read_dcdheader
734 : * Output: 0 on success, negative error code on failure.
735 : * Side effects: One timestep will be skipped; fd will be positioned at the
736 : * next timestep.
737 : */
738 0 : static int skip_dcdstep(fio_fd fd, int natoms, int nfixed, int charmm) {
739 : ptrdiff_t seekoffset = 0;
740 : ptrdiff_t rec_scale;
741 :
742 0 : if (charmm & DCD_HAS_64BIT_REC) {
743 : rec_scale=RECSCALE64BIT;
744 : } else {
745 : rec_scale=RECSCALE32BIT;
746 : }
747 :
748 : /* Skip charmm extra block */
749 0 : if ((charmm & DCD_IS_CHARMM) && (charmm & DCD_HAS_EXTRA_BLOCK)) {
750 0 : seekoffset += 4L*rec_scale + 48L + 4L*rec_scale;
751 : }
752 :
753 : /* For each atom set, seek past an int, the free atoms, and another int. */
754 0 : seekoffset += 3L * (2L*rec_scale + natoms - nfixed) * 4L;
755 :
756 : /* Assume that charmm 4th dim is the same size as the other three. */
757 0 : if ((charmm & DCD_IS_CHARMM) && (charmm & DCD_HAS_4DIMS)) {
758 0 : seekoffset += (2L*rec_scale + natoms - nfixed) * 4L;
759 : }
760 :
761 0 : if (fio_fseek(fd, seekoffset, FIO_SEEK_CUR)) return DCD_BADEOF;
762 :
763 : return DCD_SUCCESS;
764 : }
765 :
766 :
767 : /*
768 : * Write a timestep to a dcd file
769 : * Input: fd - a file struct for which a dcd header has already been written
770 : * curframe: Count of frames written to this file, starting with 1.
771 : * curstep: Count of timesteps elapsed = istart + curframe * nsavc.
772 : * natoms: number of elements in x, y, z arrays
773 : * x, y, z: pointers to atom coordinates
774 : * Output: 0 on success, negative error code on failure.
775 : * Side effects: coordinates are written to the dcd file.
776 : */
777 0 : static int write_dcdstep(fio_fd fd, int curframe, int curstep, int N,
778 : const float *X, const float *Y, const float *Z,
779 : const double *unitcell, int charmm) {
780 : int out_integer;
781 :
782 0 : if (charmm) {
783 : /* write out optional unit cell */
784 0 : if (unitcell != NULL) {
785 : out_integer = 48; /* 48 bytes (6 floats) */
786 0 : fio_write_int32(fd, out_integer);
787 0 : WRITE(fd, unitcell, out_integer);
788 0 : fio_write_int32(fd, out_integer);
789 : }
790 : }
791 :
792 : /* write out coordinates */
793 0 : out_integer = N*4; /* N*4 bytes per X/Y/Z array (N floats per array) */
794 0 : fio_write_int32(fd, out_integer);
795 0 : if (fio_fwrite((void *) X, out_integer, 1, fd) != 1) return DCD_BADWRITE;
796 0 : fio_write_int32(fd, out_integer);
797 0 : fio_write_int32(fd, out_integer);
798 0 : if (fio_fwrite((void *) Y, out_integer, 1, fd) != 1) return DCD_BADWRITE;
799 0 : fio_write_int32(fd, out_integer);
800 0 : fio_write_int32(fd, out_integer);
801 0 : if (fio_fwrite((void *) Z, out_integer, 1, fd) != 1) return DCD_BADWRITE;
802 0 : fio_write_int32(fd, out_integer);
803 :
804 : /* update the DCD header information */
805 0 : fio_fseek(fd, NFILE_POS, FIO_SEEK_SET);
806 0 : fio_write_int32(fd, curframe);
807 0 : fio_fseek(fd, NSTEP_POS, FIO_SEEK_SET);
808 0 : fio_write_int32(fd, curstep);
809 0 : fio_fseek(fd, 0, FIO_SEEK_END);
810 :
811 0 : return DCD_SUCCESS;
812 : }
813 :
814 :
815 : /*
816 : * Write a header for a new dcd file
817 : * Input: fd - file struct opened for binary writing
818 : * remarks - string to be put in the remarks section of the header.
819 : * The string will be truncated to 70 characters.
820 : * natoms, istart, nsavc, delta - see comments in read_dcdheader
821 : * Output: 0 on success, negative error code on failure.
822 : * Side effects: Header information is written to the dcd file.
823 : */
824 0 : static int write_dcdheader(fio_fd fd, const char *remarks, int N,
825 : int ISTART, int NSAVC, double DELTA, int with_unitcell,
826 : int charmm) {
827 : int out_integer;
828 : float out_float;
829 : char title_string[200];
830 : time_t cur_time;
831 : struct tm *tmbuf;
832 : char time_str[81];
833 :
834 0 : out_integer = 84;
835 0 : WRITE(fd, (char *) & out_integer, sizeof(int));
836 : strcpy(title_string, "CORD");
837 0 : WRITE(fd, title_string, 4);
838 0 : fio_write_int32(fd, 0); /* Number of frames in file, none written yet */
839 0 : fio_write_int32(fd, ISTART); /* Starting timestep */
840 0 : fio_write_int32(fd, NSAVC); /* Timesteps between frames written to the file */
841 0 : fio_write_int32(fd, 0); /* Number of timesteps in simulation */
842 0 : fio_write_int32(fd, 0); /* NAMD writes NSTEP or ISTART - NSAVC here? */
843 0 : fio_write_int32(fd, 0);
844 0 : fio_write_int32(fd, 0);
845 0 : fio_write_int32(fd, 0);
846 0 : fio_write_int32(fd, 0);
847 0 : if (charmm) {
848 0 : out_float = DELTA;
849 0 : WRITE(fd, (char *) &out_float, sizeof(float));
850 0 : if (with_unitcell) {
851 0 : fio_write_int32(fd, 1);
852 : } else {
853 0 : fio_write_int32(fd, 0);
854 : }
855 : } else {
856 0 : WRITE(fd, (char *) &DELTA, sizeof(double));
857 : }
858 0 : fio_write_int32(fd, 0);
859 0 : fio_write_int32(fd, 0);
860 0 : fio_write_int32(fd, 0);
861 0 : fio_write_int32(fd, 0);
862 0 : fio_write_int32(fd, 0);
863 0 : fio_write_int32(fd, 0);
864 0 : fio_write_int32(fd, 0);
865 0 : fio_write_int32(fd, 0);
866 0 : if (charmm) {
867 0 : fio_write_int32(fd, 24); /* Pretend to be CHARMM version 24 */
868 : } else {
869 0 : fio_write_int32(fd, 0);
870 : }
871 0 : fio_write_int32(fd, 84);
872 0 : fio_write_int32(fd, 164);
873 0 : fio_write_int32(fd, 2);
874 :
875 : strncpy(title_string, remarks, 80);
876 0 : title_string[79] = '\0';
877 0 : WRITE(fd, title_string, 80);
878 :
879 0 : cur_time=time(NULL);
880 0 : tmbuf=localtime(&cur_time);
881 0 : strftime(time_str, 80, "REMARKS Created %d %B, %Y at %R", tmbuf);
882 0 : WRITE(fd, time_str, 80);
883 :
884 0 : fio_write_int32(fd, 164);
885 0 : fio_write_int32(fd, 4);
886 0 : fio_write_int32(fd, N);
887 0 : fio_write_int32(fd, 4);
888 :
889 0 : return DCD_SUCCESS;
890 : }
891 :
892 :
893 : /*
894 : * clean up dcd data
895 : * Input: nfixed, freeind - elements as returned by read_dcdheader
896 : * Output: None
897 : * Side effects: Space pointed to by freeind is freed if necessary.
898 : */
899 1 : static void close_dcd_read(int *indexes, float *fixedcoords) {
900 1 : free(indexes);
901 1 : free(fixedcoords);
902 1 : }
903 :
904 :
905 1 : static void *open_dcd_read(const char *path, const char *filetype,
906 : int *natoms) {
907 : dcdhandle *dcd;
908 : fio_fd fd;
909 : int rc;
910 : struct stat stbuf;
911 :
912 1 : if (!path) return NULL;
913 :
914 : #if !(defined(_MSC_VER) && defined(FASTIO_NATIVEWIN32))
915 : /* See if the file exists, and get its size */
916 : memset(&stbuf, 0, sizeof(struct stat));
917 1 : if (stat(path, &stbuf)) {
918 : printf("dcdplugin) Could not access file '%s'.\n", path);
919 0 : return NULL;
920 : }
921 : #endif
922 :
923 1 : if (fio_open(path, FIO_READ, &fd) < 0) {
924 : printf("dcdplugin) Could not open file '%s' for reading.\n", path);
925 0 : return NULL;
926 : }
927 :
928 1 : dcd = (dcdhandle *)malloc(sizeof(dcdhandle));
929 : memset(dcd, 0, sizeof(dcdhandle));
930 1 : dcd->fd = fd;
931 :
932 1 : if ((rc = read_dcdheader(dcd->fd, &dcd->natoms, &dcd->nsets, &dcd->istart,
933 : &dcd->nsavc, &dcd->delta, &dcd->nfixed, &dcd->freeind,
934 : &dcd->fixedcoords, &dcd->reverse, &dcd->charmm))) {
935 0 : print_dcderror("read_dcdheader", rc);
936 0 : fio_fclose(dcd->fd);
937 0 : free(dcd);
938 0 : return NULL;
939 : }
940 :
941 : /*
942 : * Check that the file is big enough to really hold the number of sets
943 : * it claims to have. Then we'll use nsets to keep track of where EOF
944 : * should be.
945 : */
946 : {
947 : fio_size_t ndims, firstframesize, framesize, extrablocksize;
948 : fio_size_t trjsize, filesize, curpos;
949 : int newnsets;
950 :
951 1 : extrablocksize = dcd->charmm & DCD_HAS_EXTRA_BLOCK ? 48 + 8 : 0;
952 1 : ndims = dcd->charmm & DCD_HAS_4DIMS ? 4 : 3;
953 1 : firstframesize = (dcd->natoms+2) * ndims * sizeof(float) + extrablocksize;
954 1 : framesize = (dcd->natoms-dcd->nfixed+2) * ndims * sizeof(float)
955 1 : + extrablocksize;
956 :
957 : /*
958 : * It's safe to use ftell, even though ftell returns a long, because the
959 : * header size is < 4GB.
960 : */
961 :
962 1 : curpos = fio_ftell(dcd->fd); /* save current offset (end of header) */
963 :
964 : #if defined(_MSC_VER) && defined(FASTIO_NATIVEWIN32)
965 : /* the stat() call is not 64-bit savvy on Windows */
966 : /* so we have to use the fastio fseek/ftell routines for this */
967 : /* until we add a portable filesize routine for this purpose */
968 : fio_fseek(dcd->fd, 0, FIO_SEEK_END); /* seek to end of file */
969 : filesize = fio_ftell(dcd->fd);
970 : fio_fseek(dcd->fd, curpos, FIO_SEEK_SET); /* return to end of header */
971 : #else
972 1 : filesize = stbuf.st_size; /* this works ok on Unix machines */
973 : #endif
974 1 : trjsize = filesize - curpos - firstframesize;
975 1 : if (trjsize < 0) {
976 : printf("dcdplugin) file '%s' appears to contain no timesteps.\n", path);
977 0 : fio_fclose(dcd->fd);
978 0 : free(dcd);
979 0 : return NULL;
980 : }
981 :
982 1 : newnsets = trjsize / framesize + 1;
983 :
984 1 : if (dcd->nsets > 0 && newnsets != dcd->nsets) {
985 : printf("dcdplugin) Warning: DCD header claims %d frames, but \n"
986 : "dcdplugin) file size (%ld) indicates there are actually \n"
987 : "%d frames of size (%ld)\n",
988 : dcd->nsets, trjsize, newnsets, framesize);
989 : }
990 :
991 1 : dcd->nsets = newnsets;
992 1 : dcd->setsread = 0;
993 : }
994 :
995 1 : dcd->first = 1;
996 1 : dcd->x = (float *)malloc(dcd->natoms * sizeof(float));
997 1 : dcd->y = (float *)malloc(dcd->natoms * sizeof(float));
998 1 : dcd->z = (float *)malloc(dcd->natoms * sizeof(float));
999 1 : if (!dcd->x || !dcd->y || !dcd->z) {
1000 : printf("dcdplugin) Unable to allocate space for %d atoms.\n", dcd->natoms);
1001 0 : if (dcd->x)
1002 0 : free(dcd->x);
1003 0 : if (dcd->y)
1004 0 : free(dcd->y);
1005 0 : if (dcd->z)
1006 0 : free(dcd->z);
1007 0 : fio_fclose(dcd->fd);
1008 0 : free(dcd);
1009 0 : return NULL;
1010 : }
1011 1 : *natoms = dcd->natoms;
1012 1 : return dcd;
1013 : }
1014 :
1015 :
1016 22 : static int read_next_timestep(void *v, int natoms, molfile_timestep_t *ts) {
1017 : dcdhandle *dcd;
1018 : int i, j, rc;
1019 : float unitcell[6];
1020 22 : unitcell[0] = unitcell[2] = unitcell[5] = 0.0f;
1021 22 : unitcell[1] = unitcell[3] = unitcell[4] = 90.0f;
1022 : dcd = (dcdhandle *)v;
1023 :
1024 : /* Check for EOF here; that way all EOF's encountered later must be errors */
1025 22 : if (dcd->setsread == dcd->nsets) return MOLFILE_EOF;
1026 21 : dcd->setsread++;
1027 21 : if (!ts) {
1028 0 : if (dcd->first && dcd->nfixed) {
1029 : /* We can't just skip it because we need the fixed atom coordinates */
1030 0 : rc = read_dcdstep(dcd->fd, dcd->natoms, dcd->x, dcd->y, dcd->z,
1031 : unitcell, dcd->nfixed, dcd->first, dcd->freeind, dcd->fixedcoords,
1032 : dcd->reverse, dcd->charmm);
1033 0 : dcd->first = 0;
1034 0 : return rc; /* XXX this needs to be updated */
1035 : }
1036 0 : dcd->first = 0;
1037 : /* XXX this needs to be changed */
1038 0 : return skip_dcdstep(dcd->fd, dcd->natoms, dcd->nfixed, dcd->charmm);
1039 : }
1040 21 : rc = read_dcdstep(dcd->fd, dcd->natoms, dcd->x, dcd->y, dcd->z, unitcell,
1041 : dcd->nfixed, dcd->first, dcd->freeind, dcd->fixedcoords,
1042 : dcd->reverse, dcd->charmm);
1043 21 : dcd->first = 0;
1044 21 : if (rc < 0) {
1045 0 : print_dcderror("read_dcdstep", rc);
1046 0 : return MOLFILE_ERROR;
1047 : }
1048 :
1049 : /* copy timestep data from plugin-local buffers to VMD's buffer */
1050 : /* XXX
1051 : * This code is still the root of all evil. Just doing this extra copy
1052 : * cuts the I/O rate of the DCD reader from 728 MB/sec down to
1053 : * 394 MB/sec when reading from a ram filesystem.
1054 : * For a physical disk filesystem, the I/O rate goes from
1055 : * 187 MB/sec down to 122 MB/sec. Clearly this extra copy has to go.
1056 : */
1057 : {
1058 21 : int natoms = dcd->natoms;
1059 21 : float *nts = ts->coords;
1060 21 : const float *bufx = dcd->x;
1061 21 : const float *bufy = dcd->y;
1062 21 : const float *bufz = dcd->z;
1063 :
1064 483 : for (i=0, j=0; i<natoms; i++, j+=3) {
1065 462 : nts[j ] = bufx[i];
1066 462 : nts[j + 1] = bufy[i];
1067 462 : nts[j + 2] = bufz[i];
1068 : }
1069 : }
1070 :
1071 21 : ts->A = unitcell[0];
1072 21 : ts->B = unitcell[2];
1073 21 : ts->C = unitcell[5];
1074 :
1075 21 : if (unitcell[1] >= -1.0 && unitcell[1] <= 1.0 &&
1076 21 : unitcell[3] >= -1.0 && unitcell[3] <= 1.0 &&
1077 21 : unitcell[4] >= -1.0 && unitcell[4] <= 1.0) {
1078 : /* This file was generated by CHARMM, or by NAMD > 2.5, with the angle */
1079 : /* cosines of the periodic cell angles written to the DCD file. */
1080 : /* This formulation improves rounding behavior for orthogonal cells */
1081 : /* so that the angles end up at precisely 90 degrees, unlike acos(). */
1082 21 : ts->alpha = 90.0 - asin(unitcell[4]) * 90.0 / M_PI_2; /* cosBC */
1083 21 : ts->beta = 90.0 - asin(unitcell[3]) * 90.0 / M_PI_2; /* cosAC */
1084 21 : ts->gamma = 90.0 - asin(unitcell[1]) * 90.0 / M_PI_2; /* cosAB */
1085 : } else {
1086 : /* This file was likely generated by NAMD 2.5 and the periodic cell */
1087 : /* angles are specified in degrees rather than angle cosines. */
1088 0 : ts->alpha = unitcell[4]; /* angle between B and C */
1089 0 : ts->beta = unitcell[3]; /* angle between A and C */
1090 0 : ts->gamma = unitcell[1]; /* angle between A and B */
1091 : }
1092 :
1093 : return MOLFILE_SUCCESS;
1094 : }
1095 :
1096 :
1097 1 : static void close_file_read(void *v) {
1098 : dcdhandle *dcd = (dcdhandle *)v;
1099 1 : close_dcd_read(dcd->freeind, dcd->fixedcoords);
1100 1 : fio_fclose(dcd->fd);
1101 1 : free(dcd->x);
1102 1 : free(dcd->y);
1103 1 : free(dcd->z);
1104 1 : free(dcd);
1105 1 : }
1106 :
1107 :
1108 0 : static void *open_dcd_write(const char *path, const char *filetype,
1109 : int natoms) {
1110 : dcdhandle *dcd;
1111 : fio_fd fd;
1112 : int rc;
1113 : int istart, nsavc;
1114 : double delta;
1115 : int with_unitcell;
1116 : int charmm;
1117 :
1118 0 : if (fio_open(path, FIO_WRITE, &fd) < 0) {
1119 : printf("dcdplugin) Could not open file '%s' for writing\n", path);
1120 0 : return NULL;
1121 : }
1122 :
1123 0 : dcd = (dcdhandle *)malloc(sizeof(dcdhandle));
1124 : memset(dcd, 0, sizeof(dcdhandle));
1125 0 : dcd->fd = fd;
1126 :
1127 : istart = 0; /* starting timestep of DCD file */
1128 : nsavc = 1; /* number of timesteps between written DCD frames */
1129 : delta = 1.0; /* length of a timestep */
1130 :
1131 0 : if (getenv("VMDDCDWRITEXPLORFORMAT") != NULL) {
1132 : with_unitcell = 0; /* no unit cell info */
1133 : charmm = DCD_IS_XPLOR; /* X-PLOR format */
1134 : printf("dcdplugin) WARNING: Writing DCD file in X-PLOR format, \n");
1135 : printf("dcdplugin) WARNING: unit cell information will be lost!\n");
1136 : } else {
1137 : with_unitcell = 1; /* contains unit cell infor (Charmm format) */
1138 : charmm = DCD_IS_CHARMM; /* charmm-formatted DCD file */
1139 : if (with_unitcell)
1140 : charmm |= DCD_HAS_EXTRA_BLOCK;
1141 : }
1142 :
1143 0 : rc = write_dcdheader(dcd->fd, "Created by DCD plugin", natoms,
1144 : istart, nsavc, delta, with_unitcell, charmm);
1145 :
1146 0 : if (rc < 0) {
1147 0 : print_dcderror("write_dcdheader", rc);
1148 0 : fio_fclose(dcd->fd);
1149 0 : free(dcd);
1150 0 : return NULL;
1151 : }
1152 :
1153 0 : dcd->natoms = natoms;
1154 0 : dcd->nsets = 0;
1155 0 : dcd->istart = istart;
1156 0 : dcd->nsavc = nsavc;
1157 0 : dcd->with_unitcell = with_unitcell;
1158 0 : dcd->charmm = charmm;
1159 0 : dcd->x = (float *)malloc(natoms * sizeof(float));
1160 0 : dcd->y = (float *)malloc(natoms * sizeof(float));
1161 0 : dcd->z = (float *)malloc(natoms * sizeof(float));
1162 0 : return dcd;
1163 : }
1164 :
1165 :
1166 0 : static int write_timestep(void *v, const molfile_timestep_t *ts) {
1167 : dcdhandle *dcd = (dcdhandle *)v;
1168 : int i, rc, curstep;
1169 0 : float *pos = ts->coords;
1170 : double unitcell[6];
1171 0 : unitcell[0] = unitcell[2] = unitcell[5] = 0.0f;
1172 0 : unitcell[1] = unitcell[3] = unitcell[4] = 90.0f;
1173 :
1174 : /* copy atom coords into separate X/Y/Z arrays for writing */
1175 0 : for (i=0; i<dcd->natoms; i++) {
1176 0 : dcd->x[i] = *(pos++);
1177 0 : dcd->y[i] = *(pos++);
1178 0 : dcd->z[i] = *(pos++);
1179 : }
1180 0 : dcd->nsets++;
1181 0 : curstep = dcd->istart + dcd->nsets * dcd->nsavc;
1182 :
1183 0 : unitcell[0] = ts->A;
1184 0 : unitcell[2] = ts->B;
1185 0 : unitcell[5] = ts->C;
1186 0 : unitcell[1] = sin((M_PI_2 / 90.0) * (90.0 - ts->gamma)); /* cosAB */
1187 0 : unitcell[3] = sin((M_PI_2 / 90.0) * (90.0 - ts->beta)); /* cosAC */
1188 0 : unitcell[4] = sin((M_PI_2 / 90.0) * (90.0 - ts->alpha)); /* cosBC */
1189 :
1190 0 : rc = write_dcdstep(dcd->fd, dcd->nsets, curstep, dcd->natoms,
1191 0 : dcd->x, dcd->y, dcd->z,
1192 0 : dcd->with_unitcell ? unitcell : NULL,
1193 : dcd->charmm);
1194 0 : if (rc < 0) {
1195 0 : print_dcderror("write_dcdstep", rc);
1196 0 : return MOLFILE_ERROR;
1197 : }
1198 :
1199 : return MOLFILE_SUCCESS;
1200 : }
1201 :
1202 0 : static void close_file_write(void *v) {
1203 : dcdhandle *dcd = (dcdhandle *)v;
1204 0 : fio_fclose(dcd->fd);
1205 0 : free(dcd->x);
1206 0 : free(dcd->y);
1207 0 : free(dcd->z);
1208 0 : free(dcd);
1209 0 : }
1210 :
1211 :
1212 : /*
1213 : * Initialization stuff here
1214 : */
1215 : static molfile_plugin_t plugin;
1216 :
1217 10632 : VMDPLUGIN_API int VMDPLUGIN_init() {
1218 : memset(&plugin, 0, sizeof(molfile_plugin_t));
1219 10632 : plugin.abiversion = vmdplugin_ABIVERSION;
1220 10632 : plugin.type = MOLFILE_PLUGIN_TYPE;
1221 10632 : plugin.name = "dcd";
1222 10632 : plugin.prettyname = "CHARMM,NAMD,XPLOR DCD Trajectory";
1223 10632 : plugin.author = "Axel Kohlmeyer, Justin Gullingsrud, John Stone";
1224 10632 : plugin.majorv = 1;
1225 10632 : plugin.minorv = 18;
1226 10632 : plugin.is_reentrant = VMDPLUGIN_THREADSAFE;
1227 10632 : plugin.filename_extension = "dcd";
1228 10632 : plugin.open_file_read = open_dcd_read;
1229 10632 : plugin.read_next_timestep = read_next_timestep;
1230 10632 : plugin.close_file_read = close_file_read;
1231 10632 : plugin.open_file_write = open_dcd_write;
1232 10632 : plugin.write_timestep = write_timestep;
1233 10632 : plugin.close_file_write = close_file_write;
1234 10632 : return VMDPLUGIN_SUCCESS;
1235 : }
1236 :
1237 10632 : VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
1238 10632 : (*cb)(v, (vmdplugin_t *)&plugin);
1239 10632 : return VMDPLUGIN_SUCCESS;
1240 : }
1241 :
1242 0 : VMDPLUGIN_API int VMDPLUGIN_fini() {
1243 0 : return VMDPLUGIN_SUCCESS;
1244 : }
1245 :
1246 : }
1247 : }
1248 :
1249 :
1250 : #ifdef TEST_DCDPLUGIN
1251 :
1252 : #include <sys/time.h>
1253 :
1254 : /* get the time of day from the system clock, and store it (in seconds) */
1255 : double time_of_day(void) {
1256 : #if defined(_MSC_VER)
1257 : double t;
1258 :
1259 : t = GetTickCount();
1260 : t = t / 1000.0;
1261 :
1262 : return t;
1263 : #else
1264 : struct timeval tm;
1265 : struct timezone tz;
1266 :
1267 : gettimeofday(&tm, &tz);
1268 : return((double)(tm.tv_sec) + (double)(tm.tv_usec)/1000000.0);
1269 : #endif
1270 : }
1271 :
1272 : int main(int argc, char *argv[]) {
1273 : molfile_timestep_t timestep;
1274 : void *v;
1275 : dcdhandle *dcd;
1276 : int i, natoms;
1277 : float sizeMB =0.0, totalMB = 0.0;
1278 : double starttime, endtime, totaltime = 0.0;
1279 :
1280 : while (--argc) {
1281 : ++argv;
1282 : natoms = 0;
1283 : v = open_dcd_read(*argv, "dcd", &natoms);
1284 : if (!v) {
1285 : fprintf(stderr, "main) open_dcd_read failed for file %s\n", *argv);
1286 : return 1;
1287 : }
1288 : dcd = (dcdhandle *)v;
1289 : sizeMB = ((natoms * 3.0) * dcd->nsets * 4.0) / (1024.0 * 1024.0);
1290 : totalMB += sizeMB;
1291 : printf("main) file: %s\n", *argv);
1292 : printf(" %d atoms, %d frames, size: %6.1fMB\n", natoms, dcd->nsets, sizeMB);
1293 :
1294 : starttime = time_of_day();
1295 : timestep.coords = (float *)malloc(3*sizeof(float)*natoms);
1296 : for (i=0; i<dcd->nsets; i++) {
1297 : int rc = read_next_timestep(v, natoms, ×tep);
1298 : if (rc) {
1299 : fprintf(stderr, "error in read_next_timestep on frame %d\n", i);
1300 : return 1;
1301 : }
1302 : }
1303 : endtime = time_of_day();
1304 : close_file_read(v);
1305 : totaltime += endtime - starttime;
1306 : printf(" Time: %5.1f seconds\n", endtime - starttime);
1307 : printf(" Speed: %5.1f MB/sec, %5.1f timesteps/sec\n", sizeMB / (endtime - starttime), (dcd->nsets / (endtime - starttime)));
1308 : }
1309 : printf("Overall Size: %6.1f MB\n", totalMB);
1310 : printf("Overall Time: %6.1f seconds\n", totaltime);
1311 : printf("Overall Speed: %5.1f MB/sec\n", totalMB / totaltime);
1312 : return 0;
1313 : }
1314 :
1315 : #endif
1316 :
1317 : #endif
|
