1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: t; c-basic-offset: 4 -*-
2 *
3 * Copyright (c) 2009-2014, Erik Lindahl & David van der Spoel
4 * All rights reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without
7 * modification, are permitted provided that the following conditions are met:
8 *
9 * 1. Redistributions of source code must retain the above copyright notice,
10 * this
11 * list of conditions and the following disclaimer.
12 *
13 * 2. Redistributions in binary form must reproduce the above copyright notice,
14 * this list of conditions and the following disclaimer in the documentation
15 * and/or other materials provided with the distribution.
16 *
17 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 * ARE
21 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
22 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
24 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
25 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
26 * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
27 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28 */
29
30 #include "xdrfile.h"
31 #include "xdrfile_xtc.h"
32 #include <stdio.h>
33 #include <stdlib.h>
34
35 #define MAGIC 1995
36
37 enum { FALSE, TRUE };
38
xtc_header(XDRFILE * xd,int * natoms,int * step,float * time,mybool bRead)39 int xtc_header(XDRFILE *xd, int *natoms, int *step, float *time,
40 mybool bRead) {
41 int result, magic, n = 1;
42
43 /* Note: read is same as write. He he he */
44 magic = MAGIC;
45 if ((result = xdrfile_write_int(&magic, n, xd)) != n) {
46 if (bRead)
47 return exdrENDOFFILE;
48 else
49 return exdrINT;
50 }
51 if (magic != MAGIC)
52 return exdrMAGIC;
53 if ((result = xdrfile_write_int(natoms, n, xd)) != n)
54 return exdrINT;
55 if ((result = xdrfile_write_int(step, n, xd)) != n)
56 return exdrINT;
57 if ((result = xdrfile_write_float(time, n, xd)) != n)
58 return exdrFLOAT;
59
60 return exdrOK;
61 }
62
xtc_coord(XDRFILE * xd,int * natoms,matrix box,rvec * x,float * prec,mybool bRead)63 static int xtc_coord(XDRFILE *xd, int *natoms, matrix box, rvec *x, float *prec,
64 mybool bRead) {
65 int result;
66
67 /* box */
68 result = xdrfile_read_float(box[0], DIM * DIM, xd);
69 if (DIM * DIM != result)
70 return exdrFLOAT;
71 else {
72 if (bRead) {
73 result = xdrfile_decompress_coord_float(x[0], natoms, prec, xd);
74 if (result != *natoms)
75 return exdr3DX;
76 } else {
77 result = xdrfile_compress_coord_float(x[0], *natoms, *prec, xd);
78 if (result != *natoms)
79 return exdr3DX;
80 }
81 }
82 return exdrOK;
83 }
84
read_xtc_natoms(char * fn,int * natoms)85 int read_xtc_natoms(char *fn, int *natoms) {
86 XDRFILE *xd;
87 int step, result;
88 float time;
89
90 xd = xdrfile_open(fn, "r");
91 if (NULL == xd)
92 return exdrFILENOTFOUND;
93 result = xtc_header(xd, natoms, &step, &time, TRUE);
94 xdrfile_close(xd);
95
96 return result;
97 }
98
read_xtc(XDRFILE * xd,int natoms,int * step,float * time,matrix box,rvec * x,float * prec)99 int read_xtc(XDRFILE *xd, int natoms, int *step, float *time, matrix box,
100 rvec *x, float *prec)
101 /* Read subsequent frames */
102 {
103 int result;
104
105 if ((result = xtc_header(xd, &natoms, step, time, TRUE)) != exdrOK)
106 return result;
107
108 if ((result = xtc_coord(xd, &natoms, box, x, prec, 1)) != exdrOK)
109 return result;
110
111 return exdrOK;
112 }
113
114
write_xtc(XDRFILE * xd,int natoms,int step,float time,matrix box,rvec * x,float prec)115 int write_xtc(XDRFILE *xd, int natoms, int step, float time, matrix box,
116 rvec *x, float prec)
117 /* Write a frame to xtc file */
118 {
119 int result;
120
121 if ((result = xtc_header(xd, &natoms, &step, &time, FALSE)) != exdrOK)
122 return result;
123
124 if ((result = xtc_coord(xd, &natoms, box, x, &prec, 0)) != exdrOK)
125 return result;
126
127 return exdrOK;
128 }
129