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