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 : #ifndef __PLUMED_molfile_Gromacs_h
38 : #define __PLUMED_molfile_Gromacs_h
39 : /***************************************************************************
40 : *cr
41 : *cr (C) Copyright 1995-2016 The Board of Trustees of the
42 : *cr University of Illinois
43 : *cr All Rights Reserved
44 : *cr
45 : ***************************************************************************/
46 : /***************************************************************************
47 : * RCS INFORMATION:
48 : * $RCSfile: Gromacs.h,v $
49 : * $Author: johns $ $Locker: $ $State: Exp $
50 : * $Revision: 1.34 $ $Date: 2016/11/28 05:01:53 $
51 : ***************************************************************************/
52 :
53 : /*
54 : * GROMACS file format reader for VMD
55 : *
56 : * This code provides a high level I/O library for reading
57 : * and writing the following file formats:
58 : * gro GROMACS format or trajectory
59 : * g96 GROMOS-96 format or trajectory
60 : * trj Trajectory - x, v and f (binary, full precision)
61 : * trr Trajectory - x, v and f (binary, full precision, portable)
62 : * xtc Trajectory - x only (compressed, portable, any precision)
63 : * top
64 : * Currently supported: gro trj trr g96 [xtc]
65 : *
66 : * TODO list
67 : * o velocities are ignored because VMD doesn't use them, but some other
68 : * program might ...
69 : * o gro_rec() assumes positions in .gro files are nanometers and
70 : * converts to angstroms, whereas they really could be any unit
71 : */
72 :
73 : #ifndef GROMACS_H
74 : #define GROMACS_H
75 :
76 : #include <math.h>
77 : #include <stdio.h>
78 : #include <stdlib.h>
79 : #include <string.h>
80 : #include <ctype.h>
81 :
82 : #if defined(_AIX)
83 : #include <strings.h>
84 : #endif
85 :
86 : #include "endianswap.h"
87 :
88 : namespace PLMD{
89 : namespace molfile{
90 :
91 : #if defined(WIN32) || defined(WIN64)
92 : #define strcasecmp stricmp
93 : #endif
94 :
95 : #ifndef M_PI_2
96 : #define M_PI_2 1.57079632679489661922
97 : #endif
98 :
99 : // Error codes for mdio_errno
100 : #define MDIO_SUCCESS 0
101 : #define MDIO_BADFORMAT 1
102 : #define MDIO_EOF 2
103 : #define MDIO_BADPARAMS 3
104 : #define MDIO_IOERROR 4
105 : #define MDIO_BADPRECISION 5
106 : #define MDIO_BADMALLOC 6
107 : #define MDIO_CANTOPEN 7
108 : #define MDIO_BADEXTENSION 8
109 : #define MDIO_UNKNOWNFMT 9
110 : #define MDIO_CANTCLOSE 10
111 : #define MDIO_WRONGFORMAT 11
112 : #define MDIO_SIZEERROR 12
113 : #define MDIO_UNKNOWNERROR 1000
114 :
115 : #define MDIO_READ 0
116 : #define MDIO_WRITE 1
117 :
118 : #define MDIO_MAX_ERRVAL 11
119 :
120 : // Format extensions
121 : const char *mdio_fmtexts[] = {
122 : "",
123 : ".gro",
124 : ".trr",
125 : ".g96",
126 : ".trj",
127 : ".xtc",
128 : NULL
129 : };
130 :
131 :
132 : static int mdio_errcode; // Last error code
133 :
134 : #define TRX_MAGIC 1993 // Magic number for .trX files
135 : #define XTC_MAGIC 1995 // Magic number for .xtc files
136 : #define MAX_GRO_LINE 500 // Maximum line length of .gro files
137 : #define MAX_G96_LINE 500 // Maximum line length of .g96 files
138 : #define MAX_TRX_TITLE 80 // Maximum length of a title in .trX
139 : #define MAX_MDIO_TITLE 80 // Maximum supported title length
140 : #define ANGS_PER_NM 10 // Unit conversion factor
141 : #define ANGS2_PER_NM2 100 // Unit conversion factor
142 :
143 :
144 : // All the supported file types and their respective extensions
145 : #define MDFMT_GRO 1
146 : #define MDFMT_TRR 2
147 : #define MDFMT_G96 3
148 : #define MDFMT_TRJ 4
149 : #define MDFMT_XTC 5
150 :
151 :
152 : // A structure to hold .trX file format header information. This
153 : // is an optional member of the md_file structure that is used
154 : // when .trX files are being dealt with.
155 : typedef struct {
156 : int version; // File version number
157 : char title[MAX_TRX_TITLE + 1]; // File title
158 : int ir_size;
159 : int e_size;
160 : int box_size;
161 : int vir_size;
162 : int pres_size;
163 : int top_size;
164 : int sym_size;
165 : int x_size; // Positions of atoms
166 : int v_size; // Velocities of atoms
167 : int f_size;
168 : int natoms; // Number of atoms in the system
169 : int step;
170 : int nre;
171 : float t;
172 : float lambda;
173 : } trx_hdr;
174 :
175 :
176 : // A generic i/o structure that contains information about the
177 : // file itself and the input/output state
178 : typedef struct {
179 : FILE * f; // Pointer to the file
180 : int fmt; // The file format
181 : int prec; // Real number precision
182 : int rev; // Reverse endiannism?
183 : trx_hdr * trx; // Trx files require a great deal more
184 : // header data to be stored.
185 : } md_file;
186 :
187 :
188 : // A format-independent structure to hold header data from files
189 : typedef struct {
190 : char title[MAX_MDIO_TITLE + 1];
191 : int natoms;
192 : float timeval;
193 : } md_header;
194 :
195 :
196 : // A format-independent structure to hold unit cell data
197 : typedef struct {
198 : float A, B, C, alpha, beta, gamma;
199 : } md_box;
200 :
201 :
202 : // Timestep information
203 : typedef struct {
204 : float *pos; // Position array (3 * natoms)
205 : //float *vel; // Velocity array ** (VMD doesn't use this) **
206 : //float *f; // Force array ** (VMD doesn't use this) **
207 : //float *box; // Computational box ** (VMD doesn't use this) **
208 : int natoms; // Number of atoms
209 : int step; // Simulation step
210 : float time; // Time of simulation
211 : md_box *box;
212 : } md_ts;
213 :
214 :
215 : // Atom information
216 : typedef struct {
217 : char resid[7]; // Residue index number
218 : char resname[7]; // Residue name
219 : int atomnum; // Atom index number
220 : char atomname[7]; // Atom name
221 : float pos[3]; // Position array (3 * natoms)
222 : //float vel[3]; // Velocity array ** (VMD doesn't use this) **
223 : } md_atom;
224 :
225 :
226 : // Open a molecular dynamics file. The second parameter specifies
227 : // the format of the file. If it is zero, the format is determined
228 : // from the file extension. the third argument (if given) decides
229 : // whether to read (==0) or to write (!= 0).
230 : // using a default argument set to read for backward compatibility.
231 : static md_file *mdio_open(const char *, const int, const int=MDIO_READ);
232 :
233 : // Closes a molecular dynamics file.
234 : static int mdio_close(md_file *);
235 :
236 :
237 : // Format-independent file I/O routines
238 : static int mdio_header(md_file *, md_header *);
239 : static int mdio_timestep(md_file *, md_ts *);
240 :
241 :
242 : // .gro file functions
243 : static int gro_header(md_file *, char *, int, float *, int *, int = 1);
244 : static int gro_rec(md_file *, md_atom *);
245 : static int gro_timestep(md_file *, md_ts *);
246 :
247 :
248 : // .trX file functions
249 : static int trx_header(md_file *, int = 0);
250 : static int trx_int(md_file *, int *);
251 : static int trx_real(md_file *, float *);
252 :
253 : static int trx_rvector(md_file *, float *);
254 : static int trx_string(md_file *, char *, int);
255 : static int trx_timestep(md_file *, md_ts *);
256 :
257 : // .g96 file functions
258 : static int g96_header(md_file *, char *, int, float *);
259 : static int g96_timestep(md_file *, md_ts *);
260 : static int g96_rec(md_file *, md_atom *);
261 : static int g96_countatoms(md_file *);
262 :
263 :
264 : // .xtc file functions
265 : static int xtc_int(md_file *, int *);
266 : static int xtc_float(md_file *, float *);
267 : /*
268 : static int xtc_receivebits(int *, int);
269 : static void xtc_receiveints(int *, int, int, const unsigned *, int *);
270 : */
271 : static int xtc_timestep(md_file *, md_ts *);
272 : static int xtc_3dfcoord(md_file *, float *, int *, float *);
273 :
274 :
275 : // Error reporting functions
276 : static int mdio_errno(void);
277 : static const char *mdio_errmsg(int);
278 : static int mdio_seterror(int);
279 :
280 :
281 : // Miscellaneous functions
282 : static int strip_white(char *);
283 : static int mdio_readline(md_file *, char *, int, int = 1);
284 : static int mdio_tsfree(md_ts *, int = 0);
285 : static int mdio_readbox(md_box *, float *, float *, float *);
286 :
287 :
288 :
289 : static int xtc_receivebits(int *, int);
290 :
291 : // Error descriptions for mdio_errno
292 : static const char *mdio_errdescs[] = {
293 : "no error",
294 : "file does not match format",
295 : "unexpected end-of-file reached",
296 : "function called with bad parameters",
297 : "file i/o error",
298 : "unsupported precision",
299 : "out of memory",
300 : "cannot open file",
301 : "bad file extension",
302 : "unknown file format",
303 : "cannot close file",
304 : "wrong file format for this function",
305 : "binary i/o error: sizeof(int) != 4",
306 : NULL
307 : };
308 :
309 : /*! \fn static inline bool host_is_little_endian(void)
310 : * detect endiannes of host machine. returns true on little endian machines. */
311 : static inline int host_is_little_endian(void)
312 : {
313 : const union { unsigned char c[4]; unsigned int i; }
314 : fixed = { { 0x10 , 0x20 , 0x40 , 0x80 } };
315 : const unsigned int i = 0x80402010U;
316 :
317 : if (fixed.i == i) {
318 : return 1;
319 : }
320 : return 0;
321 : }
322 :
323 :
324 :
325 : // Open a molecular dynamics file. The second parameter specifies
326 : // the format of the file. If it is zero, the format is determined
327 : // from the file extension.
328 108 : md_file *mdio_open(const char *fn, const int fmt, const int rw) {
329 : md_file *mf;
330 :
331 108 : if (!fn) {
332 0 : mdio_seterror(MDIO_BADPARAMS);
333 0 : return NULL;
334 : }
335 :
336 : // Allocate memory
337 108 : mf = (md_file *) malloc(sizeof(md_file));
338 108 : if (!mf) {
339 0 : mdio_seterror(MDIO_BADMALLOC);
340 0 : return NULL;
341 : }
342 :
343 : // Zero out the structure
344 : memset(mf, 0, sizeof(md_file));
345 :
346 : // Determine the file type from the extension
347 108 : if (!fmt) {
348 : char *p;
349 : int n;
350 :
351 : // Seek to the extension part of the filename
352 0 : for (p = (char *) &fn[strlen(fn) - 1]; *p != '.' && p > fn; p--);
353 0 : if (p == fn) {
354 0 : free(mf);
355 0 : mdio_seterror(MDIO_BADEXTENSION);
356 0 : return NULL;
357 : }
358 :
359 : // Check the extension against known extensions
360 0 : for (n = 1; mdio_fmtexts[n]; n++)
361 0 : if (!strcasecmp(p, mdio_fmtexts[n])) break;
362 :
363 : // If !mdio_fmtexts[n], we failed (unknown ext)
364 0 : if (!mdio_fmtexts[n]) {
365 0 : free(mf);
366 0 : mdio_seterror(MDIO_UNKNOWNFMT);
367 0 : return NULL;
368 : }
369 :
370 : // All set
371 0 : mf->fmt = n;
372 : }
373 : else {
374 108 : mf->fmt = fmt;
375 : }
376 :
377 : // Differentiate between binary and ascii files. Also,
378 : // .trX files need a header information structure allocated.
379 108 : switch (mf->fmt) {
380 0 : case MDFMT_GRO:
381 : case MDFMT_G96: /* fallthrough */
382 0 : if (rw)
383 0 : mf->f = fopen(fn, "wt");
384 : else
385 0 : mf->f = fopen(fn, "rt");
386 :
387 : break;
388 5 : case MDFMT_TRR:
389 : case MDFMT_TRJ: /* fallthrough */
390 : // Allocate the trx header data struct
391 5 : mf->trx = (trx_hdr *) malloc(sizeof(trx_hdr));
392 5 : if (!mf->trx) {
393 0 : free(mf);
394 0 : mdio_seterror(MDIO_BADMALLOC);
395 0 : return NULL;
396 : }
397 : memset(mf->trx, 0, sizeof(trx_hdr));
398 108 : case MDFMT_XTC: /* fallthrough */
399 : // Finally, open the file
400 108 : if (rw)
401 0 : mf->f = fopen(fn, "wb");
402 : else
403 108 : mf->f = fopen(fn, "rb");
404 :
405 : break;
406 0 : default:
407 0 : free(mf);
408 0 : mdio_seterror(MDIO_UNKNOWNFMT);
409 0 : return NULL;
410 : }
411 :
412 : // Check for opening error
413 108 : if (!mf->f) {
414 0 : if (mf->trx) free(mf->trx);
415 0 : free(mf);
416 0 : mdio_seterror(MDIO_CANTOPEN);
417 0 : return NULL;
418 : }
419 :
420 : // File is opened, we're all set!
421 108 : mdio_seterror(MDIO_SUCCESS);
422 108 : return mf;
423 : }
424 :
425 :
426 : // Closes a molecular dynamics file.
427 108 : static int mdio_close(md_file *mf) {
428 108 : if (!mf) return mdio_seterror(MDIO_BADPARAMS);
429 :
430 108 : if (fclose(mf->f) == EOF) return mdio_seterror(MDIO_CANTCLOSE);
431 :
432 : // Free the dynamically allocated memory
433 108 : if (mf->trx) free(mf->trx);
434 108 : free(mf);
435 108 : return mdio_seterror(MDIO_SUCCESS);
436 : }
437 :
438 :
439 : // Returns the last error code reported by any of the mdio functions
440 108 : static int mdio_errno(void) {
441 108 : return mdio_errcode;
442 : }
443 :
444 :
445 : // Returns a textual message regarding an mdio error code
446 0 : static const char *mdio_errmsg(int n) {
447 0 : if (n < 0 || n > MDIO_MAX_ERRVAL) return (char *) "unknown error";
448 0 : else return mdio_errdescs[n];
449 : }
450 :
451 :
452 : // Sets the error code and returns an appropriate return value
453 : // for the calling function to return to its parent
454 2203476 : static int mdio_seterror(int code) {
455 2203476 : mdio_errcode = code;
456 2203476 : return code ? -1 : 0;
457 : }
458 :
459 :
460 : // Reads a line from the text file, strips leading/trailing whitespace
461 : // and newline, checks for errors, and returns the number of characters
462 : // in the string on success or -1 on error.
463 0 : static int mdio_readline(md_file *mf, char *buf, int n, int strip) {
464 0 : if (!buf || n < 1 || !mf) return mdio_seterror(MDIO_BADPARAMS);
465 :
466 : // Read the line
467 0 : fgets(buf, n, mf->f);
468 :
469 : // End of file reached?
470 0 : if (feof(mf->f)) return mdio_seterror(MDIO_EOF);
471 :
472 : // File I/O error?
473 0 : if (ferror(mf->f)) return mdio_seterror(MDIO_IOERROR);
474 :
475 : // comment line?
476 0 : if (buf[0] == '#') return mdio_readline(mf,buf,n,strip);
477 :
478 : // Strip whitespace
479 0 : if (strip) strip_white(buf);
480 :
481 0 : return strlen(buf);
482 : }
483 :
484 :
485 : // Strips leading and trailing whitespace from a string. Tabs,
486 : // spaces, newlines and carriage returns are stripped. Example:
487 : // "\n hello\t \r" becomes "hello".
488 0 : static int strip_white(char *buf) {
489 : int i, j, k;
490 :
491 : // Protect against NULL pointer
492 0 : if (!buf) return -1;
493 0 : if (!strlen(buf)) return -1;
494 :
495 : // Kill trailing whitespace first
496 0 : for (i = strlen(buf) - 1;
497 0 : buf[i] == ' ' || buf[i] == '\t' ||
498 0 : buf[i] == '\n' || buf[i] == '\r';
499 : i--)
500 0 : buf[i] = 0;
501 :
502 : // Skip past leading whitespace
503 0 : for (i = 0; buf[i] == ' ' || buf[i] == '\t' ||
504 0 : buf[i] == '\n' || buf[i] == '\r'; i++);
505 0 : if (i) {
506 : k = 0;
507 0 : for (j = i; buf[j]; j++)
508 0 : buf[k++] = buf[j];
509 0 : buf[k] = 0;
510 : }
511 :
512 0 : return strlen(buf);
513 : }
514 :
515 :
516 : // Frees the memory allocated in a ts structure. The holderror
517 : // parameter defaults to zero. Programs that are calling this
518 : // function because of an error reported by another function should
519 : // set holderror so that mdio_tsfree() does not overwrite the error
520 : // code with mdio_seterror().
521 9962 : static int mdio_tsfree(md_ts *ts, int holderror) {
522 9962 : if (!ts) {
523 0 : if (holderror) return -1;
524 0 : else return mdio_seterror(MDIO_BADPARAMS);
525 : }
526 :
527 9962 : if (ts->pos && ts->natoms > 0) free(ts->pos);
528 :
529 9962 : if (ts->box) free(ts->box);
530 :
531 9962 : if (holderror) return -1;
532 9962 : else return mdio_seterror(MDIO_SUCCESS);
533 : }
534 :
535 :
536 : // Converts box basis vectors to A, B, C, alpha, beta, and gamma.
537 : // Stores values in md_box struct, which should be allocated before calling
538 : // this function.
539 9962 : static int mdio_readbox(md_box *box, float *x, float *y, float *z) {
540 : float A, B, C;
541 :
542 9962 : if (!box) {
543 0 : return mdio_seterror(MDIO_BADPARAMS);
544 : }
545 :
546 : // A, B, C are the lengths of the x, y, z vectors, respectively
547 9962 : A = sqrt( x[0]*x[0] + x[1]*x[1] + x[2]*x[2] ) * ANGS_PER_NM;
548 9962 : B = sqrt( y[0]*y[0] + y[1]*y[1] + y[2]*y[2] ) * ANGS_PER_NM;
549 9962 : C = sqrt( z[0]*z[0] + z[1]*z[1] + z[2]*z[2] ) * ANGS_PER_NM;
550 9962 : if ((A<=0) || (B<=0) || (C<=0)) {
551 : /* Use zero-length box size and set angles to 90. */
552 0 : box->A = box->B = box->C = 0;
553 0 : box->alpha = box->beta = box->gamma = 90;
554 : } else {
555 9962 : box->A = A;
556 9962 : box->B = B;
557 9962 : box->C = C;
558 :
559 : // gamma, beta, alpha are the angles between the x & y, x & z, y & z
560 : // vectors, respectively
561 9962 : box->gamma = acos( (x[0]*y[0]+x[1]*y[1]+x[2]*y[2])*ANGS2_PER_NM2/(A*B) ) * 90.0/M_PI_2;
562 9962 : box->beta = acos( (x[0]*z[0]+x[1]*z[1]+x[2]*z[2])*ANGS2_PER_NM2/(A*C) ) * 90.0/M_PI_2;
563 9962 : box->alpha = acos( (y[0]*z[0]+y[1]*z[1]+y[2]*z[2])*ANGS2_PER_NM2/(B*C) ) * 90.0/M_PI_2;
564 : }
565 9962 : return mdio_seterror(MDIO_SUCCESS);
566 : }
567 :
568 :
569 : // Reads the header of a file (format independent)
570 108 : static int mdio_header(md_file *mf, md_header *mdh) {
571 : int n;
572 108 : if (!mf || !mdh) return mdio_seterror(MDIO_BADPARAMS);
573 108 : if (!mf->f) return mdio_seterror(MDIO_BADPARAMS);
574 :
575 108 : switch (mf->fmt) {
576 0 : case MDFMT_GRO:
577 0 : if (gro_header(mf, mdh->title, MAX_MDIO_TITLE,
578 : &mdh->timeval, &mdh->natoms, 1) < 0)
579 0 : return -1;
580 : return 0;
581 :
582 5 : case MDFMT_TRR:
583 : case MDFMT_TRJ: /* fallthrough */
584 5 : if (trx_header(mf, 1) < 0) return -1;
585 5 : mdh->natoms = mf->trx->natoms;
586 5 : mdh->timeval = (float) mf->trx->t;
587 5 : strncpy(mdh->title, mf->trx->title, MAX_MDIO_TITLE);
588 5 : return 0;
589 :
590 0 : case MDFMT_G96:
591 0 : if (g96_header(mf, mdh->title, MAX_MDIO_TITLE,
592 : &mdh->timeval) < 0) return -1;
593 0 : mdh->natoms = -1;
594 0 : return 0;
595 :
596 : case MDFMT_XTC:
597 : memset(mdh, 0, sizeof(md_header));
598 : // Check magic number
599 103 : if (xtc_int(mf, &n) < 0) return -1;
600 103 : if (n != XTC_MAGIC) return mdio_seterror(MDIO_BADFORMAT);
601 :
602 : // Get number of atoms
603 103 : if (xtc_int(mf, &n) < 0) return -1;
604 103 : mdh->natoms = n;
605 103 : rewind(mf->f);
606 : return 0;
607 :
608 0 : default:
609 0 : return mdio_seterror(MDIO_UNKNOWNFMT);
610 : }
611 : }
612 :
613 :
614 : // Reads in a timestep from a file (format independent)
615 10070 : static int mdio_timestep(md_file *mf, md_ts *ts) {
616 10070 : if (!mf || !ts) return mdio_seterror(MDIO_BADPARAMS);
617 10070 : if (!mf->f) return mdio_seterror(MDIO_BADPARAMS);
618 :
619 10070 : switch (mf->fmt) {
620 0 : case MDFMT_GRO:
621 0 : return gro_timestep(mf, ts);
622 :
623 155 : case MDFMT_TRR:
624 : case MDFMT_TRJ: /* fallthrough */
625 155 : return trx_timestep(mf, ts);
626 :
627 0 : case MDFMT_G96:
628 0 : return g96_timestep(mf, ts);
629 :
630 9915 : case MDFMT_XTC:
631 9915 : return xtc_timestep(mf, ts);
632 :
633 0 : default:
634 0 : return mdio_seterror(MDIO_UNKNOWNFMT);
635 : }
636 : }
637 :
638 :
639 :
640 0 : static int g96_header(md_file *mf, char *title, int titlelen, float *timeval) {
641 : char buf[MAX_G96_LINE + 1];
642 : char *p;
643 :
644 : // Check parameters
645 0 : if (!mf) return mdio_seterror(MDIO_BADPARAMS);
646 :
647 : // The header consists of blocks. The title block
648 : // is mandatory, and a TIMESTEP block is optional.
649 : // Example:
650 : //
651 : // TITLE
652 : // Generated by trjconv : t= 90.00000
653 : // more title info
654 : // .
655 : // .
656 : // .
657 : // END
658 : // .
659 : // .
660 : // .
661 :
662 0 : if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
663 0 : if (strcasecmp(buf, "TITLE")) return mdio_seterror(MDIO_BADFORMAT);
664 :
665 : // Read in the title itself
666 0 : if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
667 :
668 : // The timevalue can be included in the title string
669 : // after a "t=" prefix.
670 0 : if ((p = (char *) strstr(buf, "t="))) {
671 : char *q = p;
672 0 : *(q--) = 0;
673 :
674 : // Skip the `t=' and strip whitespace from
675 : // the resulting strings
676 0 : p += 2;
677 0 : strip_white(p);
678 0 : strip_white(buf);
679 :
680 : // Grab the timevalue from the title string
681 0 : if (timeval) *timeval = (float) atof(p);
682 : }
683 : else {
684 : // No timevalue - just copy the string and strip
685 : // any leading/trailing whitespace
686 0 : if (timeval) *timeval = 0;
687 0 : strip_white(buf);
688 : }
689 :
690 : // Copy the title string
691 0 : if (title && titlelen) strncpy(title, buf, titlelen);
692 :
693 : // Now ignore subsequent title lines and get the END string
694 0 : while (strcasecmp(buf, "END"))
695 0 : if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
696 :
697 : // Done!
698 0 : return mdio_seterror(MDIO_SUCCESS);
699 : }
700 :
701 :
702 : // Used to determine the number of atoms in a g96 file, because
703 : // VMD needs to know this for some reason.
704 0 : static int g96_countatoms(md_file *mf) {
705 : char buf[MAX_G96_LINE + 1];
706 : int natoms;
707 : int n;
708 : long fpos;
709 : float lastf;
710 :
711 0 : if (!mf) return mdio_seterror(MDIO_BADPARAMS);
712 :
713 0 : fpos = ftell(mf->f);
714 :
715 : natoms = 0;
716 : for (;;) {
717 0 : if (mdio_readline(mf, buf, MAX_G96_LINE + 1, 0) < 0)
718 : break;
719 0 : n = sscanf(buf, "%*6c%*6c%*6c%*6c %*f %*f %f", &lastf);
720 0 : if (n == 1) natoms++;
721 : else {
722 0 : strip_white(buf);
723 0 : if (!strcasecmp(buf, "END")) break;
724 : }
725 : }
726 :
727 0 : fseek(mf->f, fpos, SEEK_SET);
728 0 : return natoms;
729 : }
730 :
731 :
732 : // Reads an atom line from the G96 file
733 0 : static int g96_rec(md_file *mf, md_atom *ma) {
734 : char buf[MAX_G96_LINE + 1];
735 : char atomnum[7];
736 : int n;
737 :
738 : // Check parameters
739 0 : if (!mf || !ma) return mdio_seterror(MDIO_BADPARAMS);
740 :
741 : // Read in a line, assuming it is an atom line
742 : do {
743 0 : if (mdio_readline(mf, buf, MAX_G96_LINE + 1, 0) < 0) return -1;
744 0 : } while (buf[0] == '#' || strlen(buf) == 0);
745 :
746 0 : n = sscanf(buf, "%6c%6c%6c%6c %f %f %f",
747 0 : ma->resid, ma->resname, ma->atomname, atomnum,
748 : &ma->pos[0], &ma->pos[1], &ma->pos[2]);
749 0 : if (n == 7) {
750 0 : atomnum[6] = 0;
751 0 : ma->resid[6] = 0;
752 0 : ma->resname[6] = 0;
753 0 : ma->atomname[6] = 0;
754 :
755 0 : strip_white(atomnum);
756 0 : strip_white(ma->resid);
757 0 : strip_white(ma->resname);
758 0 : strip_white(ma->atomname);
759 :
760 0 : ma->atomnum = atoi(atomnum);
761 :
762 0 : ma->pos[0] *= ANGS_PER_NM;
763 0 : ma->pos[1] *= ANGS_PER_NM;
764 0 : ma->pos[2] *= ANGS_PER_NM;
765 :
766 0 : return 0;
767 : }
768 :
769 0 : return mdio_seterror(MDIO_BADFORMAT);
770 : }
771 :
772 :
773 : // Reads a timestep from a G96 file and stores the data in
774 : // the generic md_ts structure. Returns 0 on success or a
775 : // negative number on error and sets mdio_errcode.
776 0 : static int g96_timestep(md_file *mf, md_ts *ts) {
777 : char buf[MAX_G96_LINE + 1];
778 : char stripbuf[MAX_G96_LINE + 1];
779 : float pos[3], x[3], y[3], z[3], *currAtom;
780 : long fpos;
781 : int n, i, boxItems;
782 :
783 : // Check parameters
784 0 : if (!mf || !ts) return mdio_seterror(MDIO_BADPARAMS);
785 :
786 : // Allocate data space for the timestep, using the number of atoms
787 : // determined by open_g96_read().
788 0 : ts->pos = (float *) malloc(sizeof(float) * 3 * ts->natoms);
789 0 : if (!ts->pos) {
790 0 : return mdio_seterror(MDIO_BADMALLOC);
791 : }
792 : currAtom = ts->pos;
793 :
794 : // The timesteps follow the header in a fixed block
795 : // format:
796 : //
797 : // TIMESTEP
798 : // <step number> <time value>
799 : // END
800 : // POSITIONRED
801 : // <x float> <y float> <z float>
802 : // . . .
803 : // . . .
804 : // . . .
805 : // END
806 : // VELOCITYRED
807 : // <x float> <y float> <z float>
808 : // . . .
809 : // . . .
810 : // . . .
811 : // END
812 : // BOX
813 : // <x float> <y float> <z float>
814 : // END
815 : //
816 : // -----
817 : //
818 : // The TIMESTEP, VELOCITY and BOX blocks are optional.
819 : // Floats are written in 15.9 precision.
820 : //
821 : // Reference: GROMACS 2.0 user manual
822 : // http://rugmd4.chem.rug.nl/~gmx/online2.0/g96.html
823 :
824 : // First, look for an (optional) title block and skip it
825 0 : if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
826 :
827 0 : if (!strcasecmp(buf, "TITLE")) {
828 : // skip over the text until we reach 'END'
829 0 : while (strcasecmp(buf, "END")) {
830 0 : if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
831 : }
832 :
833 : // Read in the next line
834 0 : if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
835 : }
836 :
837 : // Next, look for a timestep block
838 0 : if (!strcasecmp(buf, "TIMESTEP")) {
839 : // Read in the value line
840 0 : if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
841 :
842 : // Extract the time value and the timestep index
843 0 : n = sscanf(buf, "%d %f", &ts->step, &ts->time);
844 0 : if (n != 2) return mdio_seterror(MDIO_BADFORMAT);
845 :
846 : // Read the "END" line
847 0 : if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
848 0 : if (strcasecmp(buf, "END"))
849 0 : return mdio_seterror(MDIO_BADFORMAT);
850 :
851 : // Read in the next line
852 0 : if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
853 : }
854 : else {
855 : // No timestep specified -- set to zero
856 0 : ts->step = 0;
857 0 : ts->time = 0;
858 : }
859 :
860 : // At this point a POSITION or POSITIONRED block
861 : // is REQUIRED by the format
862 0 : if (!strcasecmp(buf, "POSITIONRED")) {
863 :
864 : // So now we read in some atoms
865 : i = 0;
866 0 : while (i < ts->natoms) {
867 : // Read in an atom
868 0 : if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0)
869 : return -1;
870 :
871 : // We shouldn't reach the end yet
872 0 : if (!strcasecmp(buf, "END"))
873 0 : return mdio_seterror(MDIO_BADFORMAT);
874 :
875 : // Get the x,y,z coordinates
876 0 : n = sscanf(buf, "%f %f %f", &pos[0], &pos[1], &pos[2]);
877 :
878 : // Ignore improperly formatted lines
879 0 : if (n == 3) {
880 0 : pos[0] *= ANGS_PER_NM;
881 0 : pos[1] *= ANGS_PER_NM;
882 0 : pos[2] *= ANGS_PER_NM;
883 :
884 : // Copy the atom data into the array
885 : memcpy(currAtom, pos, sizeof(float) * 3);
886 0 : currAtom += 3;
887 0 : i++;
888 : }
889 : }
890 : }
891 0 : else if (!strcasecmp(buf, "POSITION") || !strcasecmp(buf, "REFPOSITION")) {
892 : /*
893 : char resnum[7];
894 : char resname[7];
895 : char atomname[7];
896 : char atomnum[7];
897 : */
898 :
899 : // So now we read in some atoms
900 : i = 0;
901 0 : while (i < ts->natoms) {
902 : // Read in the first line
903 0 : if (mdio_readline(mf, buf, MAX_G96_LINE + 1, 0) < 0)
904 : return -1;
905 :
906 : // We shouldn't reach the end yet
907 : strcpy(stripbuf, buf);
908 0 : strip_white(stripbuf);
909 0 : if (!strcasecmp(stripbuf, "END"))
910 0 : return mdio_seterror(MDIO_BADFORMAT);
911 :
912 : // Get the x,y,z coordinates and name data
913 0 : n = sscanf(buf, "%*6c%*6c%*6c%*6c %f %f %f",
914 : &pos[0], &pos[1], &pos[2]);
915 :
916 : // Ignore improperly formatted lines
917 0 : if (n == 3) {
918 0 : pos[0] *= ANGS_PER_NM;
919 0 : pos[1] *= ANGS_PER_NM;
920 0 : pos[2] *= ANGS_PER_NM;
921 :
922 : // Copy the atom data into the linked list item
923 : memcpy(currAtom, pos, sizeof(float) * 3);
924 0 : currAtom += 3;
925 0 : i++;
926 : }
927 : }
928 : }
929 : else {
930 0 : return mdio_seterror(MDIO_BADFORMAT);
931 : }
932 :
933 : // Read the END keyword
934 0 : if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0)
935 : return -1;
936 0 : if (strcasecmp(buf, "END"))
937 0 : return mdio_seterror(MDIO_BADFORMAT);
938 :
939 : // ... another problem: there may or may not be a VELOCITY
940 : // block or a BOX block, so we need to read one line beyond
941 : // the POSITION block to determine this. If neither VEL. nor
942 : // BOX are present we've read a line too far and infringed
943 : // on the next timestep, so we need to keep track of the
944 : // position now for a possible fseek() later to backtrack.
945 0 : fpos = ftell(mf->f);
946 :
947 : // Now we must read in the velocities and the box, if present
948 0 : if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) {
949 : // It's okay if we end the file here; any other errors need to be
950 : // reported.
951 0 : if (mdio_errcode == MDIO_EOF)
952 0 : return mdio_seterror(MDIO_SUCCESS);
953 : else
954 : return -1;
955 : }
956 :
957 : // Is there a velocity block present ?
958 0 : if (!strcasecmp(buf, "VELOCITY") || !strcasecmp(buf, "VELOCITYRED")) {
959 : // Ignore all the coordinates - VMD doesn't use them
960 : for (;;) {
961 0 : if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0)
962 : return -1;
963 0 : if (!strcasecmp(buf, "END")) break;
964 : }
965 :
966 : // Again, record our position because we may need
967 : // to fseek here later if we read too far.
968 0 : fpos = ftell(mf->f);
969 :
970 : // Go ahead and read the next line.
971 0 : if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
972 : }
973 :
974 : // Is there a box present ?
975 0 : if (!strcasecmp(buf, "BOX")) {
976 0 : if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) return -1;
977 0 : boxItems = sscanf(buf, " %f %f %f %f %f %f %f %f %f",
978 : &x[0], &y[1], &z[2], &x[1], &x[2], &y[0], &y[2], &z[0], &z[1]);
979 0 : if (boxItems == 3) {
980 0 : x[1] = x[2] = 0;
981 0 : y[0] = y[2] = 0;
982 0 : z[0] = z[1] = 0;
983 : }
984 0 : else if (boxItems != 9)
985 0 : return mdio_seterror(MDIO_BADFORMAT);
986 :
987 : // Allocate the box and convert the vectors.
988 0 : ts->box = (md_box *) malloc(sizeof(md_box));
989 0 : if (mdio_readbox(ts->box, x, y, z) < 0) {
990 0 : free(ts->box);
991 0 : ts->box = NULL;
992 0 : return mdio_seterror(MDIO_BADFORMAT);
993 : }
994 :
995 0 : if (mdio_readline(mf, buf, MAX_G96_LINE + 1) < 0) {
996 0 : free(ts->box);
997 0 : ts->box = NULL;
998 0 : return -1;
999 : }
1000 0 : if (strcasecmp(buf, "END")) {
1001 0 : free(ts->box);
1002 0 : ts->box = NULL;
1003 0 : return mdio_seterror(MDIO_BADFORMAT);
1004 : }
1005 : }
1006 : else {
1007 : // We have read too far, so fseek back to the
1008 : // last known safe position so we don't return
1009 : // with the file pointer set infringing on the
1010 : // next timestep data.
1011 0 : fseek(mf->f, fpos, SEEK_SET);
1012 : }
1013 :
1014 : // We're done!
1015 0 : return mdio_seterror(MDIO_SUCCESS);
1016 : }
1017 :
1018 :
1019 : // Attempts to read header data from a GROMACS structure file
1020 : // The GROMACS header format is as follows (fixed, 2 lines ASCII):
1021 : // <title> [ n= <timevalue> ]
1022 : // <num atoms>
1023 0 : static int gro_header(md_file *mf, char *title, int titlelen, float *timeval,
1024 : int *natoms, int rewind) {
1025 : char buf[MAX_GRO_LINE + 1];
1026 : long fpos;
1027 : char *p;
1028 :
1029 : // Check parameters
1030 0 : if (!mf)
1031 0 : return mdio_seterror(MDIO_BADPARAMS);
1032 :
1033 : // Get the current file position for rewinding later
1034 0 : fpos = ftell(mf->f);
1035 :
1036 : // The header consists of 2 lines - get the first line
1037 0 : if (mdio_readline(mf, buf, MAX_GRO_LINE + 1) < 0) return -1;
1038 :
1039 : // The timevalue can be included in the title string
1040 : // after a "t=" prefix.
1041 0 : if ((p = (char *) strstr(buf, "t="))) {
1042 : char *q = p;
1043 0 : *(q--) = 0;
1044 :
1045 : // Skip the `t=' and strip whitespace from
1046 : // the resulting strings
1047 0 : p += 2;
1048 0 : strip_white(p);
1049 0 : strip_white(buf);
1050 :
1051 : // Grab the timevalue from the title string
1052 0 : if (timeval) *timeval = (float) atof(p);
1053 : } else {
1054 : // No timevalue - just copy the string
1055 0 : if (timeval) *timeval = 0;
1056 : }
1057 :
1058 : // Copy the title string
1059 0 : if (title && titlelen) strncpy(title, buf, titlelen);
1060 :
1061 : // Get the second line and grab the number of atoms
1062 0 : if (mdio_readline(mf, buf, MAX_GRO_LINE + 1) < 0) return -1;
1063 :
1064 : // Store the number of atoms
1065 0 : if (natoms && (!(*natoms = atoi(buf))))
1066 0 : return mdio_seterror(MDIO_BADFORMAT);
1067 :
1068 : // Now we rewind the file so that subsequent calls to
1069 : // gro_timestep() will succeed. gro_timestep() requires
1070 : // the header to be at the current file pointer.
1071 0 : if (rewind)
1072 0 : fseek(mf->f, fpos, SEEK_SET);
1073 :
1074 : return 0; // Done!
1075 : }
1076 :
1077 :
1078 : // Reads one atom record from a GROMACS file. Returns GMX_SUCCESS
1079 : // on success or a negative number on error.
1080 : //
1081 : // Record format (one line, fixed):
1082 : // rrrrrRRRRRaaaaaAAAAA <x pos> <y pos> <z pos> <x vel> <y vel> <z vel>
1083 : //
1084 : // r = residue number
1085 : // R = residue name
1086 : // a = atom name
1087 : // A = atom number
1088 : //
1089 0 : static int gro_rec(md_file *mf, md_atom *ma) {
1090 : char buf[MAX_GRO_LINE + 1];
1091 : char atomnum[6];
1092 : char xposc[12], yposc[12], zposc[12];
1093 : int n;
1094 :
1095 0 : if (!mf)
1096 0 : return mdio_seterror(MDIO_BADPARAMS);
1097 :
1098 : do {
1099 0 : if (mdio_readline(mf, buf, MAX_GRO_LINE + 1, 0) < 0)
1100 : return -1;
1101 0 : } while (buf[0] == '#' || !strlen(buf));
1102 :
1103 : // Read in the fields
1104 0 : n = sscanf(buf, "%5c%5c%5c%5c%8c%8c%8c",
1105 0 : ma->resid, ma->resname, ma->atomname, atomnum,
1106 : xposc, yposc, zposc);
1107 :
1108 0 : if (n != 7)
1109 0 : return mdio_seterror(MDIO_BADFORMAT);
1110 :
1111 : // Null terminate the strings
1112 0 : ma->resname[5] = 0;
1113 0 : ma->atomname[5] = 0;
1114 0 : ma->resid[5] = 0;
1115 0 : atomnum[5] = 0;
1116 0 : xposc[8] = 0;
1117 0 : yposc[8] = 0;
1118 0 : zposc[8] = 0;
1119 :
1120 0 : if ((sscanf(xposc, "%f", &ma->pos[0]) != 1) ||
1121 0 : (sscanf(yposc, "%f", &ma->pos[1]) != 1) ||
1122 0 : (sscanf(zposc, "%f", &ma->pos[2]) != 1)) {
1123 0 : return mdio_seterror(MDIO_BADFORMAT);
1124 : }
1125 :
1126 : // Convert strings to numbers
1127 0 : strip_white(atomnum);
1128 0 : ma->atomnum = atoi(atomnum);
1129 :
1130 : // Convert nanometers to angstroms
1131 0 : ma->pos[0] *= ANGS_PER_NM;
1132 0 : ma->pos[1] *= ANGS_PER_NM;
1133 0 : ma->pos[2] *= ANGS_PER_NM;
1134 :
1135 : // Strip leading and trailing whitespace
1136 0 : strip_white(ma->atomname);
1137 0 : strip_white(ma->resname);
1138 0 : strip_white(ma->resid);
1139 :
1140 0 : return 0;
1141 : }
1142 :
1143 :
1144 : // Reads in a timestep from a .gro file. Ignores the data
1145 : // not needed for a timestep, so is a little faster than
1146 : // calling gro_rec() for each atom. Also reads in the
1147 : // header block.
1148 : //
1149 0 : static int gro_timestep(md_file *mf, md_ts *ts) {
1150 : char buf[MAX_GRO_LINE + 1];
1151 : long coord;
1152 : int i, n, boxItems;
1153 : float x[3], y[3], z[3];
1154 : char xposc[12], yposc[12], zposc[12];
1155 :
1156 0 : if (!mf || !ts)
1157 0 : return mdio_seterror(MDIO_BADPARAMS);
1158 :
1159 0 : if (gro_header(mf, NULL, 0, &ts->time, &ts->natoms, 0) < 0)
1160 : return -1;
1161 :
1162 0 : ts->pos = (float *) malloc(3 * sizeof(float) * ts->natoms);
1163 0 : if (!ts->pos)
1164 0 : return mdio_seterror(MDIO_BADMALLOC);
1165 :
1166 : coord = 0;
1167 0 : for (i = 0; i < ts->natoms; i++) {
1168 0 : if (mdio_readline(mf, buf, MAX_GRO_LINE + 1, 0) < 0) {
1169 0 : free(ts->pos);
1170 0 : return -1;
1171 : }
1172 :
1173 0 : n = sscanf(buf, "%*5c%*5c%*5c%*5c%8c%8c%8c", xposc, yposc, zposc);
1174 0 : if (n != 3)
1175 0 : return mdio_seterror(MDIO_BADFORMAT);
1176 :
1177 0 : if ((sscanf(xposc, "%f", &ts->pos[coord ]) != 1) ||
1178 0 : (sscanf(yposc, "%f", &ts->pos[coord + 1]) != 1) ||
1179 0 : (sscanf(zposc, "%f", &ts->pos[coord + 2]) != 1)) {
1180 0 : return mdio_seterror(MDIO_BADFORMAT);
1181 : }
1182 :
1183 0 : ts->pos[coord ] *= ANGS_PER_NM;
1184 0 : ts->pos[coord + 1] *= ANGS_PER_NM;
1185 0 : ts->pos[coord + 2] *= ANGS_PER_NM;
1186 :
1187 0 : coord += 3;
1188 : }
1189 :
1190 : // Read the box, stored as three vectors representing its edges
1191 0 : if (mdio_readline(mf, buf, MAX_GRO_LINE + 1, 0) < 0) {
1192 0 : free(ts->pos);
1193 0 : return -1;
1194 : }
1195 :
1196 0 : boxItems = sscanf(buf, " %f %f %f %f %f %f %f %f %f",
1197 : &x[0], &y[1], &z[2], &x[1], &x[2], &y[0], &y[2], &z[0], &z[1]);
1198 :
1199 : // File may only include three scalars for the box information -- if
1200 : // that's the case, the box is orthoganal.
1201 0 : if (boxItems == 3) {
1202 0 : x[1] = x[2] = 0;
1203 0 : y[0] = y[2] = 0;
1204 0 : z[0] = z[1] = 0;
1205 0 : } else if (boxItems != 9) {
1206 0 : free(ts->pos);
1207 0 : return -1;
1208 : }
1209 :
1210 : // Allocate the box and convert the vectors.
1211 0 : ts->box = (md_box *) malloc(sizeof(md_box));
1212 0 : if (mdio_readbox(ts->box, x, y, z) < 0) {
1213 0 : free(ts->pos);
1214 0 : free(ts->box);
1215 0 : ts->box = NULL;
1216 0 : return -1;
1217 : }
1218 :
1219 : return 0;
1220 : }
1221 :
1222 :
1223 : // Attempts to read header data from a .trX trajectory file
1224 : //
1225 : // The .trX header format is as follows:
1226 : //
1227 : // 4 bytes - magic number (0x07C9)
1228 : // ...
1229 : //
1230 160 : static int trx_header(md_file *mf, int rewind) {
1231 : int magic;
1232 : trx_hdr *hdr;
1233 : long fpos;
1234 :
1235 160 : if (!mf) return mdio_seterror(MDIO_BADPARAMS);
1236 :
1237 : // In case we need to rewind
1238 160 : fpos = ftell(mf->f);
1239 :
1240 : // We need to store some data to the trX header data
1241 : // structure inside the md_file structure
1242 160 : hdr = mf->trx;
1243 160 : if (!mf->trx) return mdio_seterror(MDIO_BADPARAMS);
1244 :
1245 : // Read the magic number
1246 160 : if (trx_int(mf, &magic) < 0) return -1;
1247 155 : if (magic != TRX_MAGIC) {
1248 : // Try reverse endianism
1249 5 : swap4_aligned(&magic, 1);
1250 5 : if (magic != TRX_MAGIC) return mdio_seterror(MDIO_BADFORMAT);
1251 :
1252 : // Enable byte swapping (actually works, too!)
1253 5 : mf->rev = 1;
1254 : }
1255 :
1256 : // Read the version number.
1257 : // XXX. this is not the version number, but the storage size
1258 : // of the following XDR encoded string.
1259 : // the 'title' string is in fact the version identifier.
1260 : // since VMD does not use any of that, it does no harm,
1261 : // but is should still be fixed occasionally. AK 2005/01/08.
1262 :
1263 155 : if(mf->fmt!=MDFMT_TRJ) {
1264 : // It appears that TRJ files either don't contain a version
1265 : // number or don't have a length-delimiter on the string,
1266 : // whereas TRR files do contain both. Thus, with TRJ, we just
1267 : // assume that the version number is the string length and
1268 : // just hope for the best. -- WLD 2006/07/09
1269 155 : if (trx_int(mf, &hdr->version) < 0) return -1;
1270 : }
1271 :
1272 : // Read in the title string
1273 155 : if (trx_string(mf, hdr->title, MAX_TRX_TITLE) < 0)
1274 : return -1;
1275 :
1276 : // Read in some size data
1277 155 : if (trx_int(mf, &hdr->ir_size) < 0) return -1;
1278 155 : if (trx_int(mf, &hdr->e_size) < 0) return -1;
1279 155 : if (trx_int(mf, &hdr->box_size) < 0) return -1;
1280 155 : if (trx_int(mf, &hdr->vir_size) < 0) return -1;
1281 155 : if (trx_int(mf, &hdr->pres_size) < 0) return -1;
1282 155 : if (trx_int(mf, &hdr->top_size) < 0) return -1;
1283 155 : if (trx_int(mf, &hdr->sym_size) < 0) return -1;
1284 155 : if (trx_int(mf, &hdr->x_size) < 0) return -1;
1285 155 : if (trx_int(mf, &hdr->v_size) < 0) return -1;
1286 155 : if (trx_int(mf, &hdr->f_size) < 0) return -1;
1287 155 : if (trx_int(mf, &hdr->natoms) < 0) return -1;
1288 155 : if (trx_int(mf, &hdr->step) < 0) return -1;
1289 155 : if (trx_int(mf, &hdr->nre) < 0) return -1;
1290 :
1291 : // Make sure there are atoms...
1292 155 : if (!hdr->natoms) return mdio_seterror(MDIO_BADFORMAT);
1293 :
1294 : // Try to determine precision (float? double?)
1295 155 : if (hdr->x_size) mf->prec = hdr->x_size / (hdr->natoms * 3);
1296 0 : else if (hdr->v_size) mf->prec = hdr->v_size / (hdr->natoms * 3);
1297 0 : else if (hdr->f_size) mf->prec = hdr->f_size / (hdr->natoms * 3);
1298 0 : else return mdio_seterror(MDIO_BADPRECISION);
1299 :
1300 155 : if (mf->prec != sizeof(float) && mf->prec != sizeof(double)) {
1301 : // We have no data types this size! The
1302 : // file must've been generated on another
1303 : // platform
1304 0 : return mdio_seterror(MDIO_BADPRECISION);
1305 : }
1306 :
1307 : // Read in timestep and lambda
1308 155 : if (trx_real(mf, &hdr->t) < 0) return -1;
1309 155 : if (trx_real(mf, &hdr->lambda) < 0) return -1;
1310 :
1311 : // Rewind if necessary
1312 155 : if (rewind) fseek(mf->f, fpos, SEEK_SET);
1313 :
1314 : return 0;
1315 : }
1316 :
1317 :
1318 : // Reads in an integer and stores it in y. Returns GMX_SUCCESS
1319 : // on success or a negative number on error.
1320 2485 : static int trx_int(md_file *mf, int *y) {
1321 2485 : if (!mf) return mdio_seterror(MDIO_BADPARAMS);
1322 :
1323 : // sanity check.
1324 : if (sizeof(int) != 4) return mdio_seterror(MDIO_SIZEERROR);
1325 :
1326 2485 : if (y) {
1327 4970 : if (fread(y, 4, 1, mf->f) != 1)
1328 5 : return mdio_seterror(MDIO_IOERROR);
1329 2480 : if (mf->rev) swap4_aligned(y, 1);
1330 : }
1331 0 : else if (fseek(mf->f, 4, SEEK_CUR) != 0)
1332 0 : return mdio_seterror(MDIO_IOERROR);
1333 :
1334 2480 : return mdio_seterror(MDIO_SUCCESS);
1335 : }
1336 :
1337 :
1338 : // Reads in either a float or a double, depending on the
1339 : // precision, and stores that number in y. Returns
1340 : // GMX_SUCCESS on success or a negative number on error.
1341 1458760 : static int trx_real(md_file *mf, float *y) {
1342 : double x;
1343 :
1344 1458760 : if (!mf) return mdio_seterror(MDIO_BADPARAMS);
1345 :
1346 1458760 : switch (mf->prec) {
1347 1458760 : case sizeof(float):
1348 1458760 : if (!y) {
1349 0 : if (fseek(mf->f, mf->prec, SEEK_CUR) != 0)
1350 0 : return mdio_seterror(MDIO_IOERROR);
1351 : } else {
1352 2917520 : if (fread(y, mf->prec, 1, mf->f) != 1)
1353 0 : return mdio_seterror(MDIO_IOERROR);
1354 1458760 : if (mf->rev) swap4_aligned(y, 1);
1355 : }
1356 1458760 : return mdio_seterror(MDIO_SUCCESS);
1357 :
1358 0 : case sizeof(double):
1359 0 : if (!y) {
1360 0 : if (fseek(mf->f, mf->prec, SEEK_CUR) != 0)
1361 0 : return mdio_seterror(MDIO_IOERROR);
1362 : } else {
1363 0 : if (fread(&x, mf->prec, 1, mf->f) != 1)
1364 0 : return mdio_seterror(MDIO_IOERROR);
1365 0 : if (mf->rev) swap8_aligned(&x, 1);
1366 0 : *y = (float) x;
1367 : }
1368 0 : return mdio_seterror(MDIO_SUCCESS);
1369 :
1370 0 : default:
1371 0 : return mdio_seterror(MDIO_BADPRECISION);
1372 : }
1373 :
1374 : }
1375 :
1376 :
1377 : // Reads in a real-valued vector (taking precision into account).
1378 : // Stores the vector in vec, and returns GMX_SUCCESS on success
1379 : // or a negative number on error.
1380 486150 : static int trx_rvector(md_file *mf, float *vec) {
1381 486150 : if (!mf) return mdio_seterror(MDIO_BADPARAMS);
1382 :
1383 486150 : if (!vec) {
1384 0 : if (trx_real(mf, NULL) < 0) return -1;
1385 0 : if (trx_real(mf, NULL) < 0) return -1;
1386 0 : if (trx_real(mf, NULL) < 0) return -1;
1387 0 : return mdio_seterror(MDIO_SUCCESS);
1388 : } else {
1389 486150 : if (trx_real(mf, &vec[0]) < 0) return -1;
1390 486150 : if (trx_real(mf, &vec[1]) < 0) return -1;
1391 486150 : if (trx_real(mf, &vec[2]) < 0) return -1;
1392 486150 : return mdio_seterror(MDIO_SUCCESS);
1393 : }
1394 : }
1395 :
1396 :
1397 : // Reads in a string by first reading an integer containing the
1398 : // string's length, then reading in the string itself and storing
1399 : // it in str. If the length is greater than max, it is truncated
1400 : // and the rest of the string is skipped in the file. Returns the
1401 : // length of the string on success or a negative number on error.
1402 155 : static int trx_string(md_file *mf, char *str, int max) {
1403 : int size;
1404 : size_t ssize;
1405 :
1406 155 : if (!mf) return mdio_seterror(MDIO_BADPARAMS);
1407 :
1408 155 : if (trx_int(mf, &size) < 0) return -1;
1409 155 : ssize = (size_t)size;
1410 :
1411 155 : if (str && size <= max) {
1412 310 : if (fread(str, 1, size, mf->f) != ssize)
1413 0 : return mdio_seterror(MDIO_IOERROR);
1414 155 : str[size] = 0;
1415 155 : return size;
1416 0 : } else if (str) {
1417 0 : if (fread(str, 1, max, mf->f) != ssize)
1418 0 : return mdio_seterror(MDIO_IOERROR);
1419 0 : if (fseek(mf->f, size - max, SEEK_CUR) != 0)
1420 0 : return mdio_seterror(MDIO_IOERROR);
1421 0 : str[max] = 0;
1422 0 : return max;
1423 : } else {
1424 0 : if (fseek(mf->f, size, SEEK_CUR) != 0)
1425 0 : return mdio_seterror(MDIO_IOERROR);
1426 : return 0;
1427 : }
1428 : }
1429 :
1430 :
1431 : // Reads in a timestep frame from the .trX file and returns the
1432 : // data in a timestep structure. Returns NULL on error.
1433 155 : static int trx_timestep(md_file *mf, md_ts *ts) {
1434 : int i;
1435 : float x[3], y[3], z[3];
1436 : trx_hdr *hdr;
1437 :
1438 155 : if (!mf || !ts) return mdio_seterror(MDIO_BADPARAMS);
1439 155 : if (mf->fmt != MDFMT_TRJ && mf->fmt != MDFMT_TRR)
1440 0 : return mdio_seterror(MDIO_WRONGFORMAT);
1441 :
1442 : // Read the header
1443 155 : if (trx_header(mf) < 0) return -1;
1444 :
1445 : // We need some data from the trX header
1446 150 : hdr = mf->trx;
1447 150 : if (!hdr) return mdio_seterror(MDIO_BADPARAMS);
1448 :
1449 150 : if (hdr->box_size) { // XXX need to check value of box_size!!
1450 150 : if (trx_rvector(mf, x) < 0) return -1;
1451 150 : if (trx_rvector(mf, y) < 0) return -1;
1452 150 : if (trx_rvector(mf, z) < 0) return -1;
1453 : // Allocate the box and convert the vectors.
1454 150 : ts->box = (md_box *) malloc(sizeof(md_box));
1455 150 : if (mdio_readbox(ts->box, x, y, z) < 0) {
1456 0 : free(ts->box);
1457 0 : ts->box = NULL;
1458 0 : return -1;
1459 : }
1460 : }
1461 :
1462 150 : if (hdr->vir_size) {
1463 0 : if (trx_rvector(mf, NULL) < 0) return -1;
1464 0 : if (trx_rvector(mf, NULL) < 0) return -1;
1465 0 : if (trx_rvector(mf, NULL) < 0) return -1;
1466 : }
1467 :
1468 150 : if (hdr->pres_size) {
1469 0 : if (trx_rvector(mf, NULL) < 0) return -1;
1470 0 : if (trx_rvector(mf, NULL) < 0) return -1;
1471 0 : if (trx_rvector(mf, NULL) < 0) return -1;
1472 : }
1473 :
1474 150 : if (hdr->x_size) {
1475 150 : ts->pos = (float *) malloc(sizeof(float) * 3 * hdr->natoms);
1476 150 : if (!ts->pos) return mdio_seterror(MDIO_BADMALLOC);
1477 :
1478 150 : ts->natoms = hdr->natoms;
1479 :
1480 485850 : for (i = 0; i < hdr->natoms; i++) {
1481 485700 : if (trx_rvector(mf, &ts->pos[i * 3]) < 0) {
1482 0 : mdio_tsfree(ts, 1);
1483 0 : return -1;
1484 : }
1485 485700 : ts->pos[i * 3] *= ANGS_PER_NM;
1486 485700 : ts->pos[i * 3 + 1] *= ANGS_PER_NM;
1487 485700 : ts->pos[i * 3 + 2] *= ANGS_PER_NM;
1488 : }
1489 : }
1490 :
1491 150 : if (hdr->v_size) {
1492 0 : for (i = 0; i < hdr->natoms; i++) {
1493 0 : if (trx_rvector(mf, NULL) < 0) {
1494 0 : mdio_tsfree(ts, 1);
1495 0 : return -1;
1496 : }
1497 : }
1498 : }
1499 :
1500 150 : if (hdr->f_size) {
1501 0 : for (i = 0; i < hdr->natoms; i++) {
1502 0 : if (trx_rvector(mf, NULL) < 0) {
1503 0 : mdio_tsfree(ts, 1);
1504 0 : return -1;
1505 : }
1506 : }
1507 : }
1508 :
1509 150 : return mdio_seterror(MDIO_SUCCESS);
1510 : }
1511 :
1512 :
1513 : // writes an int in big endian. Returns GMX_SUCCESS
1514 : // on success or a negative number on error.
1515 0 : static int put_trx_int(md_file *mf, int y) {
1516 0 : if (!mf) return mdio_seterror(MDIO_BADPARAMS);
1517 :
1518 : // sanity check.
1519 : if (sizeof(int) != 4) return mdio_seterror(MDIO_SIZEERROR);
1520 :
1521 0 : if (mf->rev) swap4_aligned(&y, 1);
1522 0 : if (fwrite(&y, 4, 1, mf->f) != 1)
1523 0 : return mdio_seterror(MDIO_IOERROR);
1524 :
1525 0 : return mdio_seterror(MDIO_SUCCESS);
1526 : }
1527 :
1528 : // writes a real in big-endian. Returns GMX_SUCCESS
1529 : // on success or a negative number on error.
1530 0 : static int put_trx_real(md_file *mf, float y) {
1531 0 : if (!mf) return mdio_seterror(MDIO_BADPARAMS);
1532 :
1533 0 : if (mf->rev) swap4_aligned(&y, 1);
1534 0 : if (fwrite(&y, 4, 1, mf->f) != 1)
1535 0 : return mdio_seterror(MDIO_IOERROR);
1536 :
1537 0 : return mdio_seterror(MDIO_SUCCESS);
1538 : }
1539 :
1540 :
1541 : // writes an xdr encoded string. Returns GMX_SUCCESS
1542 : // on success or a negative number on error.
1543 0 : static int put_trx_string(md_file *mf, const char *s) {
1544 0 : if (!mf || !s) return mdio_seterror(MDIO_BADPARAMS);
1545 :
1546 : // write: size of object, string length, string data
1547 0 : size_t len = strlen(s);
1548 0 : if ( put_trx_int(mf, len+1)
1549 0 : || put_trx_int(mf, len)
1550 0 : || (fwrite(s, len, 1, mf->f) != 1))
1551 0 : return mdio_seterror(MDIO_IOERROR);
1552 :
1553 0 : return mdio_seterror(MDIO_SUCCESS);
1554 : }
1555 :
1556 :
1557 : // xtc_int() - reads an integer from an xtc file
1558 118045 : static int xtc_int(md_file *mf, int *i) {
1559 : unsigned char c[4];
1560 :
1561 118045 : if (!mf) return mdio_seterror(MDIO_BADPARAMS);
1562 : // sanity check.
1563 : if (sizeof(int) != 4) return mdio_seterror(MDIO_SIZEERROR);
1564 :
1565 118045 : if (fread(c, 1, 4, mf->f) != 4) {
1566 103 : if (feof(mf->f)) return mdio_seterror(MDIO_EOF);
1567 0 : else if (ferror(mf->f)) return mdio_seterror(MDIO_IOERROR);
1568 0 : else return mdio_seterror(MDIO_UNKNOWNERROR);
1569 : }
1570 :
1571 117942 : if (i) *i = c[3] + (c[2] << 8) + (c[1] << 16) + (c[0] << 24);
1572 117942 : return mdio_seterror(MDIO_SUCCESS);
1573 : }
1574 :
1575 :
1576 : // xtc_float() - reads a float from an xtc file
1577 107934 : static int xtc_float(md_file *mf, float *f) {
1578 : unsigned char c[4];
1579 : int i;
1580 :
1581 107934 : if (!mf) return mdio_seterror(MDIO_BADPARAMS);
1582 :
1583 107934 : if (fread(c, 1, 4, mf->f) != 4) {
1584 0 : if (feof(mf->f)) return mdio_seterror(MDIO_EOF);
1585 0 : else if (ferror(mf->f)) return mdio_seterror(MDIO_IOERROR);
1586 0 : else return mdio_seterror(MDIO_UNKNOWNERROR);
1587 : }
1588 :
1589 107934 : if (f) {
1590 : // By reading the number in as an integer and then
1591 : // copying it to a floating point number we can
1592 : // ensure proper endianness
1593 107934 : i = c[3] + (c[2] << 8) + (c[1] << 16) + (c[0] << 24);
1594 : memcpy(f, &i, 4);
1595 : }
1596 107934 : return mdio_seterror(MDIO_SUCCESS);
1597 : }
1598 :
1599 :
1600 : // xtc_data() - reads a specific amount of data from an xtc
1601 : // file using the xdr format.
1602 9811 : static int xtc_data(md_file *mf, char *buf, int len) {
1603 9811 : if (!mf || len < 1) return mdio_seterror(MDIO_BADPARAMS);
1604 9811 : size_t slen = (size_t)len;
1605 9811 : if (buf) {
1606 19622 : if (fread(buf, 1, slen, mf->f) != slen) {
1607 0 : if (feof(mf->f)) return mdio_seterror(MDIO_EOF);
1608 0 : if (ferror(mf->f)) return mdio_seterror(MDIO_IOERROR);
1609 0 : else return mdio_seterror(MDIO_UNKNOWNERROR);
1610 : }
1611 9811 : if (len % 4) {
1612 7617 : if (fseek(mf->f, 4 - (len % 4), SEEK_CUR)) {
1613 0 : if (feof(mf->f)) return mdio_seterror(MDIO_EOF);
1614 0 : if (ferror(mf->f)) return mdio_seterror(MDIO_IOERROR);
1615 0 : else return mdio_seterror(MDIO_UNKNOWNERROR);
1616 : }
1617 : }
1618 : }
1619 : else {
1620 : int newlen;
1621 : newlen = len;
1622 0 : if (len % 4) newlen += (4 - (len % 4));
1623 0 : if (fseek(mf->f, newlen, SEEK_CUR)) {
1624 0 : if (feof(mf->f)) return mdio_seterror(MDIO_EOF);
1625 0 : if (ferror(mf->f)) return mdio_seterror(MDIO_IOERROR);
1626 0 : else return mdio_seterror(MDIO_UNKNOWNERROR);
1627 : }
1628 : }
1629 : return len;
1630 : }
1631 :
1632 :
1633 : // xtc_timestep() - reads a timestep from an .xtc file.
1634 9915 : static int xtc_timestep(md_file *mf, md_ts *ts) {
1635 : int n;
1636 : float f, x[3], y[3], z[3];
1637 :
1638 9915 : int size = 0; // explicitly initialized to zero.
1639 : float precision;
1640 :
1641 9915 : if (!mf || !ts) return mdio_seterror(MDIO_BADPARAMS);
1642 9915 : if (!mf->f) return mdio_seterror(MDIO_BADPARAMS);
1643 9915 : if (mf->fmt != MDFMT_XTC) return mdio_seterror(MDIO_WRONGFORMAT);
1644 :
1645 : // Check magic number
1646 9915 : if (xtc_int(mf, &n) < 0) return -1;
1647 9812 : if (n != XTC_MAGIC) return mdio_seterror(MDIO_BADFORMAT);
1648 :
1649 : // Get number of atoms
1650 9812 : if (xtc_int(mf, &n) < 0) return -1;
1651 9812 : ts->natoms = n;
1652 :
1653 : // Get the simulation step
1654 9812 : if (xtc_int(mf, &n) < 0) return -1;
1655 9812 : ts->step = n;
1656 :
1657 : // Get the time value
1658 9812 : if (xtc_float(mf, &f) < 0) return -1;
1659 9812 : ts->time = f;
1660 :
1661 : // Read the basis vectors of the box
1662 19624 : if ( (xtc_float(mf, &x[0]) < 0) ||
1663 19624 : (xtc_float(mf, &x[1]) < 0) ||
1664 19624 : (xtc_float(mf, &x[2]) < 0) ||
1665 19624 : (xtc_float(mf, &y[0]) < 0) ||
1666 19624 : (xtc_float(mf, &y[1]) < 0) ||
1667 19624 : (xtc_float(mf, &y[2]) < 0) ||
1668 19624 : (xtc_float(mf, &z[0]) < 0) ||
1669 29436 : (xtc_float(mf, &z[1]) < 0) ||
1670 9812 : (xtc_float(mf, &z[2]) < 0) )
1671 0 : return -1;
1672 : // Allocate the box and convert the vectors.
1673 9812 : ts->box = (md_box *) malloc(sizeof(md_box));
1674 9812 : if (mdio_readbox(ts->box, x, y, z) < 0) {
1675 0 : free(ts->box);
1676 0 : ts->box = NULL;
1677 0 : return -1;
1678 : }
1679 :
1680 9812 : ts->pos = (float *) malloc(sizeof(float) * 3 * ts->natoms);
1681 9812 : if (!ts->pos) return mdio_seterror(MDIO_BADMALLOC);
1682 9812 : n = xtc_3dfcoord(mf, ts->pos, &size, &precision);
1683 9812 : if (n < 0) return -1;
1684 :
1685 : /* Now we're left with the job of scaling... */
1686 3959549 : for (n = 0; n < ts->natoms * 3; n++)
1687 3949737 : ts->pos[n] *= ANGS_PER_NM;
1688 :
1689 9812 : return mdio_seterror(MDIO_SUCCESS);
1690 : }
1691 :
1692 :
1693 : ///////////////////////////////////////////////////////////////////////
1694 : // This algorithm is an implementation of the 3dfcoord algorithm
1695 : // written by Frans van Hoesel (hoesel@chem.rug.nl) as part of the
1696 : // Europort project in 1995.
1697 : ///////////////////////////////////////////////////////////////////////
1698 :
1699 : // integer table used in decompression
1700 : static int xtc_magicints[] = {
1701 : 0, 0, 0, 0, 0, 0, 0, 0, 0,8, 10, 12, 16, 20, 25, 32, 40, 50, 64,
1702 : 80, 101, 128, 161, 203, 256, 322, 406, 512, 645, 812, 1024, 1290,
1703 : 1625, 2048, 2580, 3250, 4096, 5060, 6501, 8192, 10321, 13003, 16384,
1704 : 20642, 26007, 32768, 41285, 52015, 65536, 82570, 104031, 131072,
1705 : 165140, 208063, 262144, 330280, 416127, 524287, 660561, 832255,
1706 : 1048576, 1321122, 1664510, 2097152, 2642245, 3329021, 4194304,
1707 : 5284491, 6658042, 8388607, 10568983, 13316085, 16777216 };
1708 :
1709 : #define FIRSTIDX 9
1710 : /* note that magicints[FIRSTIDX-1] == 0 */
1711 : #define LASTIDX (sizeof(xtc_magicints) / sizeof(*xtc_magicints))
1712 :
1713 :
1714 : // returns the number of bits in the binary expansion of
1715 : // the given integer.
1716 0 : static int xtc_sizeofint(int size) {
1717 : unsigned int num = 1;
1718 0 : unsigned int ssize = (unsigned int)size;
1719 : int nbits = 0;
1720 :
1721 0 : while (ssize >= num && nbits < 32) {
1722 0 : nbits++;
1723 0 : num <<= 1;
1724 : }
1725 0 : return nbits;
1726 : }
1727 :
1728 : // calculates the number of bits a set of integers, when compressed,
1729 : // will take up.
1730 9811 : static int xtc_sizeofints(int nints, unsigned int *sizes) {
1731 : int i;
1732 : unsigned int num;
1733 : unsigned int nbytes, nbits, bytes[32], bytecnt, tmp;
1734 : nbytes = 1;
1735 9811 : bytes[0] = 1;
1736 : nbits = 0;
1737 39244 : for (i=0; i < nints; i++) {
1738 : tmp = 0;
1739 88504 : for (bytecnt = 0; bytecnt < nbytes; bytecnt++) {
1740 59071 : tmp = bytes[bytecnt] * sizes[i] + tmp;
1741 59071 : bytes[bytecnt] = tmp & 0xff;
1742 59071 : tmp >>= 8;
1743 : }
1744 61740 : while (tmp != 0) {
1745 32307 : bytes[bytecnt++] = tmp & 0xff;
1746 32307 : tmp >>= 8;
1747 : }
1748 : nbytes = bytecnt;
1749 : }
1750 : num = 1;
1751 9811 : nbytes--;
1752 43837 : while (bytes[nbytes] >= num) {
1753 34026 : nbits++;
1754 34026 : num *= 2;
1755 : }
1756 9811 : return nbits + nbytes * 8;
1757 : }
1758 :
1759 : // reads bits from a buffer.
1760 5721918 : static int xtc_receivebits(int *buf, int nbits) {
1761 : int cnt, num;
1762 : unsigned int lastbits, lastbyte;
1763 : unsigned char * cbuf;
1764 5721918 : int mask = (1 << nbits) -1;
1765 :
1766 5721918 : cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
1767 5721918 : cnt = buf[0];
1768 5721918 : lastbits = (unsigned int) buf[1];
1769 5721918 : lastbyte = (unsigned int) buf[2];
1770 :
1771 : num = 0;
1772 9713261 : while (nbits >= 8) {
1773 3991343 : lastbyte = ( lastbyte << 8 ) | cbuf[cnt++];
1774 3991343 : num |= (lastbyte >> lastbits) << (nbits - 8);
1775 : nbits -=8;
1776 : }
1777 5721918 : if (nbits > 0) {
1778 1730575 : if (lastbits < (unsigned int)nbits) {
1779 971927 : lastbits += 8;
1780 971927 : lastbyte = (lastbyte << 8) | cbuf[cnt++];
1781 : }
1782 1730575 : lastbits -= nbits;
1783 1730575 : num |= (lastbyte >> lastbits) & ((1 << nbits) -1);
1784 : }
1785 5721918 : num &= mask;
1786 5721918 : buf[0] = cnt;
1787 5721918 : buf[1] = lastbits;
1788 5721918 : buf[2] = lastbyte;
1789 5721918 : return num;
1790 : }
1791 :
1792 : // decompresses small integers from the buffer
1793 : // sizes parameter has to be non-zero to prevent divide-by-zero
1794 1316578 : static void xtc_receiveints(int *buf, const int nints, int nbits,
1795 : unsigned int *sizes, int *nums) {
1796 : int bytes[32];
1797 : int i, j, nbytes, p, num;
1798 :
1799 1316578 : bytes[1] = bytes[2] = bytes[3] = 0;
1800 : nbytes = 0;
1801 5268463 : while (nbits > 8) {
1802 3951885 : bytes[nbytes++] = xtc_receivebits(buf, 8);
1803 3951885 : nbits -= 8;
1804 : }
1805 1316578 : if (nbits > 0) {
1806 1316578 : bytes[nbytes++] = xtc_receivebits(buf, nbits);
1807 : }
1808 3949734 : for (i = nints-1; i > 0; i--) {
1809 : num = 0;
1810 13170082 : for (j = nbytes-1; j >=0; j--) {
1811 10536926 : num = (num << 8) | bytes[j];
1812 10536926 : p = num / sizes[i];
1813 10536926 : bytes[j] = p;
1814 10536926 : num = num - p * sizes[i];
1815 : }
1816 2633156 : nums[i] = num;
1817 : }
1818 1316578 : nums[0] = bytes[0] | (bytes[1] << 8) | (bytes[2] << 16) | (bytes[3] << 24);
1819 1316578 : }
1820 :
1821 : // function that actually reads and writes compressed coordinates
1822 9812 : static int xtc_3dfcoord(md_file *mf, float *fp, int *size, float *precision) {
1823 : static int *ip = NULL;
1824 : static int oldsize;
1825 : static int *buf;
1826 :
1827 : int minint[3], maxint[3], *lip;
1828 : int smallidx;
1829 : unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3;
1830 : int flag, k;
1831 : int small, smaller, i, is_smaller, run;
1832 : float *lfp;
1833 : int tmp, *thiscoord, prevcoord[3];
1834 :
1835 : int bufsize, lsize;
1836 : unsigned int bitsize;
1837 : float inv_precision;
1838 :
1839 : /* avoid uninitialized data compiler warnings */
1840 : bitsizeint[0] = 0;
1841 : bitsizeint[1] = 0;
1842 : bitsizeint[2] = 0;
1843 :
1844 9812 : if (xtc_int(mf, &lsize) < 0) return -1;
1845 :
1846 9812 : if (*size != 0 && lsize != *size) return mdio_seterror(MDIO_BADFORMAT);
1847 :
1848 9812 : *size = lsize;
1849 9812 : size3 = *size * 3;
1850 9812 : if (*size <= 9) {
1851 2 : for (i = 0; i < *size; i++) {
1852 1 : if (xtc_float(mf, fp + (3 * i)) < 0) return -1;
1853 1 : if (xtc_float(mf, fp + (3 * i) + 1) < 0) return -1;
1854 1 : if (xtc_float(mf, fp + (3 * i) + 2) < 0) return -1;
1855 : }
1856 : return *size;
1857 : }
1858 9811 : xtc_float(mf, precision);
1859 9811 : if (ip == NULL) {
1860 102 : ip = (int *)malloc(size3 * sizeof(*ip));
1861 102 : if (ip == NULL) return mdio_seterror(MDIO_BADMALLOC);
1862 102 : bufsize = (int) (size3 * 1.2);
1863 102 : buf = (int *)malloc(bufsize * sizeof(*buf));
1864 102 : if (buf == NULL) return mdio_seterror(MDIO_BADMALLOC);
1865 102 : oldsize = *size;
1866 9709 : } else if (*size > oldsize) {
1867 0 : ip = (int *)realloc(ip, size3 * sizeof(*ip));
1868 0 : if (ip == NULL) return mdio_seterror(MDIO_BADMALLOC);
1869 0 : bufsize = (int) (size3 * 1.2);
1870 0 : buf = (int *)realloc(buf, bufsize * sizeof(*buf));
1871 0 : if (buf == NULL) return mdio_seterror(MDIO_BADMALLOC);
1872 0 : oldsize = *size;
1873 : }
1874 9811 : buf[0] = buf[1] = buf[2] = 0;
1875 :
1876 9811 : xtc_int(mf, &(minint[0]));
1877 9811 : xtc_int(mf, &(minint[1]));
1878 9811 : xtc_int(mf, &(minint[2]));
1879 :
1880 9811 : xtc_int(mf, &(maxint[0]));
1881 9811 : xtc_int(mf, &(maxint[1]));
1882 9811 : xtc_int(mf, &(maxint[2]));
1883 :
1884 9811 : sizeint[0] = maxint[0] - minint[0]+1;
1885 9811 : sizeint[1] = maxint[1] - minint[1]+1;
1886 9811 : sizeint[2] = maxint[2] - minint[2]+1;
1887 :
1888 : /* check if one of the sizes is to big to be multiplied */
1889 9811 : if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) {
1890 0 : bitsizeint[0] = xtc_sizeofint(sizeint[0]);
1891 0 : bitsizeint[1] = xtc_sizeofint(sizeint[1]);
1892 0 : bitsizeint[2] = xtc_sizeofint(sizeint[2]);
1893 : bitsize = 0; /* flag the use of large sizes */
1894 : } else {
1895 9811 : bitsize = xtc_sizeofints(3, sizeint);
1896 : }
1897 :
1898 9811 : xtc_int(mf, &smallidx);
1899 9811 : smaller = xtc_magicints[FIRSTIDX > smallidx - 1 ? FIRSTIDX : smallidx - 1] / 2;
1900 9811 : small = xtc_magicints[smallidx] / 2;
1901 9811 : sizesmall[0] = sizesmall[1] = sizesmall[2] = xtc_magicints[smallidx];
1902 :
1903 : /* check for zero values that would yield corrupted data */
1904 9811 : if ( !sizesmall[0] || !sizesmall[1] || !sizesmall[2] ) {
1905 : printf("XTC corrupted, sizesmall==0 (case 1)\n");
1906 0 : return -1;
1907 : }
1908 :
1909 :
1910 : /* buf[0] holds the length in bytes */
1911 9811 : if (xtc_int(mf, &(buf[0])) < 0) return -1;
1912 :
1913 9811 : if (xtc_data(mf, (char *) &buf[3], (int) buf[0]) < 0) return -1;
1914 :
1915 9811 : buf[0] = buf[1] = buf[2] = 0;
1916 :
1917 : lfp = fp;
1918 9811 : inv_precision = 1.0f / (*precision);
1919 : run = 0;
1920 : i = 0;
1921 9811 : lip = ip;
1922 301460 : while (i < lsize) {
1923 291649 : thiscoord = (int *)(lip) + i * 3;
1924 :
1925 291649 : if (bitsize == 0) {
1926 0 : thiscoord[0] = xtc_receivebits(buf, bitsizeint[0]);
1927 0 : thiscoord[1] = xtc_receivebits(buf, bitsizeint[1]);
1928 0 : thiscoord[2] = xtc_receivebits(buf, bitsizeint[2]);
1929 : } else {
1930 291649 : xtc_receiveints(buf, 3, bitsize, sizeint, thiscoord);
1931 : }
1932 :
1933 291649 : i++;
1934 291649 : thiscoord[0] += minint[0];
1935 291649 : thiscoord[1] += minint[1];
1936 291649 : thiscoord[2] += minint[2];
1937 :
1938 : prevcoord[0] = thiscoord[0];
1939 : prevcoord[1] = thiscoord[1];
1940 : prevcoord[2] = thiscoord[2];
1941 :
1942 :
1943 291649 : flag = xtc_receivebits(buf, 1);
1944 : is_smaller = 0;
1945 291649 : if (flag == 1) {
1946 161806 : run = xtc_receivebits(buf, 5);
1947 161806 : is_smaller = run % 3;
1948 161806 : run -= is_smaller;
1949 161806 : is_smaller--;
1950 : }
1951 291649 : if (run > 0) {
1952 239992 : thiscoord += 3;
1953 1264921 : for (k = 0; k < run; k+=3) {
1954 1024929 : xtc_receiveints(buf, 3, smallidx, sizesmall, thiscoord);
1955 1024929 : i++;
1956 1024929 : thiscoord[0] += prevcoord[0] - small;
1957 1024929 : thiscoord[1] += prevcoord[1] - small;
1958 1024929 : thiscoord[2] += prevcoord[2] - small;
1959 1024929 : if (k == 0) {
1960 : /* interchange first with second atom for better
1961 : * compression of water molecules
1962 : */
1963 239992 : tmp = thiscoord[0]; thiscoord[0] = prevcoord[0];
1964 : prevcoord[0] = tmp;
1965 239992 : tmp = thiscoord[1]; thiscoord[1] = prevcoord[1];
1966 : prevcoord[1] = tmp;
1967 239992 : tmp = thiscoord[2]; thiscoord[2] = prevcoord[2];
1968 : prevcoord[2] = tmp;
1969 239992 : *lfp++ = prevcoord[0] * inv_precision;
1970 239992 : *lfp++ = prevcoord[1] * inv_precision;
1971 239992 : *lfp++ = prevcoord[2] * inv_precision;
1972 :
1973 239992 : if ( !sizesmall[0] || !sizesmall[1] || !sizesmall[2] ) {
1974 : printf("XTC corrupted, sizesmall==0 (case 2)\n");
1975 0 : return -1;
1976 : }
1977 :
1978 : } else {
1979 : prevcoord[0] = thiscoord[0];
1980 : prevcoord[1] = thiscoord[1];
1981 : prevcoord[2] = thiscoord[2];
1982 : }
1983 1024929 : *lfp++ = thiscoord[0] * inv_precision;
1984 1024929 : *lfp++ = thiscoord[1] * inv_precision;
1985 1024929 : *lfp++ = thiscoord[2] * inv_precision;
1986 : }
1987 : } else {
1988 51657 : *lfp++ = thiscoord[0] * inv_precision;
1989 51657 : *lfp++ = thiscoord[1] * inv_precision;
1990 51657 : *lfp++ = thiscoord[2] * inv_precision;
1991 : }
1992 291649 : smallidx += is_smaller;
1993 291649 : if (is_smaller < 0) {
1994 : small = smaller;
1995 32417 : if (smallidx > FIRSTIDX) {
1996 32417 : smaller = xtc_magicints[smallidx - 1] /2;
1997 : } else {
1998 : smaller = 0;
1999 : }
2000 259232 : } else if (is_smaller > 0) {
2001 : smaller = small;
2002 96566 : small = xtc_magicints[smallidx] / 2;
2003 : }
2004 291649 : sizesmall[0] = sizesmall[1] = sizesmall[2] = xtc_magicints[smallidx] ;
2005 : }
2006 : return 1;
2007 : }
2008 : }
2009 : }
2010 : #endif
2011 : #endif
|
