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 <stdio.h>
31 #include <stdlib.h>
32 #include <string.h>
33
34 #ifdef HAVE_CONFIG_H
35 #include "config.h"
36 #endif
37
38 #include "xdrfile.h"
39 #include "xdrfile_trr.h"
40
41 #define BUFSIZE 128
42 #define GROMACS_MAGIC 1993
43
nFloatSize(t_trnheader * sh,int * nflsz)44 static int nFloatSize(t_trnheader *sh, int *nflsz) {
45 int nflsize = 0;
46
47 if (sh->box_size)
48 nflsize = sh->box_size / (DIM * DIM);
49 else if (sh->x_size)
50 nflsize = sh->x_size / (sh->natoms * DIM);
51 else if (sh->v_size)
52 nflsize = sh->v_size / (sh->natoms * DIM);
53 else if (sh->f_size)
54 nflsize = sh->f_size / (sh->natoms * DIM);
55 else
56 return exdrHEADER;
57
58 if (((nflsize != sizeof(float)) && (nflsize != sizeof(double))))
59 return exdrHEADER;
60
61 *nflsz = nflsize;
62
63 return exdrOK;
64 }
65
do_trnheader(XDRFILE * xd,mybool bRead,t_trnheader * sh)66 int do_trnheader(XDRFILE *xd, mybool bRead, t_trnheader *sh) {
67 int magic = GROMACS_MAGIC;
68 int nflsz, slen, result;
69 char *version = "GMX_trn_file";
70 char buf[BUFSIZE];
71
72 if (xdrfile_read_int(&magic, 1, xd) != 1)
73 return exdrINT;
74
75 if (bRead) {
76 if (xdrfile_read_int(&slen, 1, xd) != 1)
77 return exdrINT;
78 if (slen != strlen(version) + 1)
79 return exdrSTRING;
80 if (xdrfile_read_string(buf, BUFSIZE, xd) <= 0)
81 return exdrSTRING;
82 } else {
83 slen = strlen(version) + 1;
84 if (xdrfile_read_int(&slen, 1, xd) != 1)
85 return exdrINT;
86 if (xdrfile_write_string(version, xd) != (strlen(version) + 1))
87 return exdrSTRING;
88 }
89 if (xdrfile_read_int(&sh->ir_size, 1, xd) != 1)
90 return exdrINT;
91 if (xdrfile_read_int(&sh->e_size, 1, xd) != 1)
92 return exdrINT;
93 if (xdrfile_read_int(&sh->box_size, 1, xd) != 1)
94 return exdrINT;
95 if (xdrfile_read_int(&sh->vir_size, 1, xd) != 1)
96 return exdrINT;
97 if (xdrfile_read_int(&sh->pres_size, 1, xd) != 1)
98 return exdrINT;
99 if (xdrfile_read_int(&sh->top_size, 1, xd) != 1)
100 return exdrINT;
101 if (xdrfile_read_int(&sh->sym_size, 1, xd) != 1)
102 return exdrINT;
103 if (xdrfile_read_int(&sh->x_size, 1, xd) != 1)
104 return exdrINT;
105 if (xdrfile_read_int(&sh->v_size, 1, xd) != 1)
106 return exdrINT;
107 if (xdrfile_read_int(&sh->f_size, 1, xd) != 1)
108 return exdrINT;
109 if (xdrfile_read_int(&sh->natoms, 1, xd) != 1)
110 return exdrINT;
111
112 if ((result = nFloatSize(sh, &nflsz)) != exdrOK)
113 return result;
114 sh->bDouble = (nflsz == sizeof(double));
115
116 if (xdrfile_read_int(&sh->step, 1, xd) != 1)
117 return exdrINT;
118 if (xdrfile_read_int(&sh->nre, 1, xd) != 1)
119 return exdrINT;
120 if (sh->bDouble) {
121 if (xdrfile_read_double(&sh->td, 1, xd) != 1)
122 return exdrDOUBLE;
123 sh->tf = sh->td;
124 if (xdrfile_read_double(&sh->lambdad, 1, xd) != 1)
125 return exdrDOUBLE;
126 sh->lambdaf = sh->lambdad;
127 } else {
128 if (xdrfile_read_float(&sh->tf, 1, xd) != 1)
129 return exdrFLOAT;
130 sh->td = sh->tf;
131 if (xdrfile_read_float(&sh->lambdaf, 1, xd) != 1)
132 return exdrFLOAT;
133 sh->lambdad = sh->lambdaf;
134 }
135
136 return exdrOK;
137 }
138
do_htrn(XDRFILE * xd,mybool bRead,t_trnheader * sh,matrix box,rvec * x,rvec * v,rvec * f)139 static int do_htrn(XDRFILE *xd, mybool bRead, t_trnheader *sh, matrix box,
140 rvec *x, rvec *v, rvec *f) {
141 double pvd[DIM * DIM];
142 double *dx = NULL;
143 float pvf[DIM * DIM];
144 float *fx = NULL;
145 int i, j;
146
147 if (sh->bDouble) {
148 if (sh->box_size != 0) {
149 if (!bRead) {
150 for (i = 0; (i < DIM); i++)
151 for (j = 0; (j < DIM); j++)
152 if (NULL != box) {
153 pvd[i * DIM + j] = box[i][j];
154 }
155 }
156 if (xdrfile_read_double(pvd, DIM * DIM, xd) == DIM * DIM) {
157 for (i = 0; (i < DIM); i++)
158 for (j = 0; (j < DIM); j++)
159 if (NULL != box) {
160 box[i][j] = pvd[i * DIM + j];
161 }
162 } else
163 return exdrDOUBLE;
164 }
165
166 if (sh->vir_size != 0) {
167 if (xdrfile_read_double(pvd, DIM * DIM, xd) != DIM * DIM)
168 return exdrDOUBLE;
169 }
170
171 if (sh->pres_size != 0) {
172 if (xdrfile_read_double(pvd, DIM * DIM, xd) != DIM * DIM)
173 return exdrDOUBLE;
174 }
175
176 if ((sh->x_size != 0) || (sh->v_size != 0) || (sh->f_size != 0)) {
177 dx = (double *)calloc(sh->natoms * DIM, sizeof(dx[0]));
178 if (NULL == dx)
179 return exdrNOMEM;
180 }
181 if (sh->x_size != 0) {
182 if (!bRead) {
183 for (i = 0; (i < sh->natoms); i++)
184 for (j = 0; (j < DIM); j++)
185 if (NULL != x) {
186 dx[i * DIM + j] = x[i][j];
187 }
188 }
189 if (xdrfile_read_double(dx, sh->natoms * DIM, xd) == sh->natoms * DIM) {
190 if (bRead) {
191 for (i = 0; (i < sh->natoms); i++)
192 for (j = 0; (j < DIM); j++)
193 if (NULL != x) {
194 x[i][j] = dx[i * DIM + j];
195 }
196 }
197 } else
198 return exdrDOUBLE;
199 }
200 if (sh->v_size != 0) {
201 if (!bRead) {
202 for (i = 0; (i < sh->natoms); i++)
203 for (j = 0; (j < DIM); j++)
204 if (NULL != x) {
205 dx[i * DIM + j] = v[i][j];
206 }
207 }
208 if (xdrfile_read_double(dx, sh->natoms * DIM, xd) == sh->natoms * DIM) {
209 for (i = 0; (i < sh->natoms); i++)
210 for (j = 0; (j < DIM); j++)
211 if (NULL != v) {
212 v[i][j] = dx[i * DIM + j];
213 }
214 } else
215 return exdrDOUBLE;
216 }
217 if (sh->f_size != 0) {
218 if (!bRead) {
219 for (i = 0; (i < sh->natoms); i++)
220 for (j = 0; (j < DIM); j++)
221 if (NULL != x) {
222 dx[i * DIM + j] = f[i][j];
223 }
224 }
225 if (xdrfile_read_double(dx, sh->natoms * DIM, xd) == sh->natoms * DIM) {
226 for (i = 0; (i < sh->natoms); i++) {
227 for (j = 0; (j < DIM); j++) {
228 if (NULL != f) {
229 f[i][j] = dx[i * DIM + j];
230 }
231 }
232 }
233 } else
234 return exdrDOUBLE;
235 }
236 if ((sh->x_size != 0) || (sh->v_size != 0) || (sh->f_size != 0)) {
237 free(dx);
238 }
239 } else
240 /* Float */
241 {
242 if (sh->box_size != 0) {
243 if (!bRead) {
244 for (i = 0; (i < DIM); i++)
245 for (j = 0; (j < DIM); j++)
246 if (NULL != box) {
247 pvf[i * DIM + j] = box[i][j];
248 }
249 }
250 if (xdrfile_read_float(pvf, DIM * DIM, xd) == DIM * DIM) {
251 for (i = 0; (i < DIM); i++) {
252 for (j = 0; (j < DIM); j++) {
253 if (NULL != box) {
254 box[i][j] = pvf[i * DIM + j];
255 }
256 }
257 }
258 } else
259 return exdrFLOAT;
260 }
261
262 if (sh->vir_size != 0) {
263 if (xdrfile_read_float(pvf, DIM * DIM, xd) != DIM * DIM)
264 return exdrFLOAT;
265 }
266
267 if (sh->pres_size != 0) {
268 if (xdrfile_read_float(pvf, DIM * DIM, xd) != DIM * DIM)
269 return exdrFLOAT;
270 }
271
272 if ((sh->x_size != 0) || (sh->v_size != 0) || (sh->f_size != 0)) {
273 fx = (float *)calloc(sh->natoms * DIM, sizeof(fx[0]));
274 if (NULL == fx)
275 return exdrNOMEM;
276 }
277 if (sh->x_size != 0) {
278 if (!bRead) {
279 for (i = 0; (i < sh->natoms); i++)
280 for (j = 0; (j < DIM); j++)
281 if (NULL != x) {
282 fx[i * DIM + j] = x[i][j];
283 }
284 }
285 if (xdrfile_read_float(fx, sh->natoms * DIM, xd) == sh->natoms * DIM) {
286 if (bRead) {
287 for (i = 0; (i < sh->natoms); i++)
288 for (j = 0; (j < DIM); j++)
289 if (NULL != x)
290 x[i][j] = fx[i * DIM + j];
291 }
292 } else
293 return exdrFLOAT;
294 }
295 if (sh->v_size != 0) {
296 if (!bRead) {
297 for (i = 0; (i < sh->natoms); i++)
298 for (j = 0; (j < DIM); j++)
299 if (NULL != x) {
300 fx[i * DIM + j] = v[i][j];
301 }
302 }
303 if (xdrfile_read_float(fx, sh->natoms * DIM, xd) == sh->natoms * DIM) {
304 for (i = 0; (i < sh->natoms); i++)
305 for (j = 0; (j < DIM); j++)
306 if (NULL != v)
307 v[i][j] = fx[i * DIM + j];
308 } else
309 return exdrFLOAT;
310 }
311 if (sh->f_size != 0) {
312 if (!bRead) {
313 for (i = 0; (i < sh->natoms); i++)
314 for (j = 0; (j < DIM); j++)
315 if (NULL != x) {
316 fx[i * DIM + j] = f[i][j];
317 }
318 }
319 if (xdrfile_read_float(fx, sh->natoms * DIM, xd) == sh->natoms * DIM) {
320 for (i = 0; (i < sh->natoms); i++)
321 for (j = 0; (j < DIM); j++)
322 if (NULL != f)
323 f[i][j] = fx[i * DIM + j];
324 } else
325 return exdrFLOAT;
326 }
327 if ((sh->x_size != 0) || (sh->v_size != 0) || (sh->f_size != 0)) {
328 free(fx);
329 }
330 }
331 return exdrOK;
332 }
333
do_trn(XDRFILE * xd,mybool bRead,int * step,float * t,float * lambda,matrix box,int * natoms,rvec * x,rvec * v,rvec * f,int * has_prop)334 static int do_trn(XDRFILE *xd, mybool bRead, int *step, float *t, float *lambda,
335 matrix box, int *natoms, rvec *x, rvec *v, rvec *f,
336 int *has_prop) {
337 t_trnheader *sh;
338 int result;
339
340 sh = (t_trnheader *)calloc(1, sizeof(*sh));
341
342 if (!bRead) {
343 sh->box_size = (NULL != box) ? sizeof(matrix) : 0;
344 sh->x_size = ((NULL != x) ? (*natoms * sizeof(x[0])) : 0);
345 sh->v_size = ((NULL != v) ? (*natoms * sizeof(v[0])) : 0);
346 sh->f_size = ((NULL != f) ? (*natoms * sizeof(f[0])) : 0);
347 sh->natoms = *natoms;
348 sh->step = *step;
349 sh->nre = 0;
350 sh->td = *t;
351 sh->lambdad = *lambda;
352 sh->tf = *t;
353 sh->lambdaf = *lambda;
354 }
355 if ((result = do_trnheader(xd, bRead, sh)) != exdrOK)
356 return result;
357 if (bRead) {
358 *natoms = sh->natoms;
359 *step = sh->step;
360 *t = sh->td;
361 *lambda = sh->lambdad;
362 /* Flag what we read */
363 *has_prop = 0;
364 if (sh->x_size > 0) {
365 *has_prop |= HASX;
366 }
367 if (sh->v_size > 0) {
368 *has_prop |= HASV;
369 }
370 if (sh->f_size > 0) {
371 *has_prop |= HASF;
372 }
373 }
374 if ((result = do_htrn(xd, bRead, sh, box, x, v, f)) != exdrOK)
375 return result;
376
377 free(sh);
378
379 return exdrOK;
380 }
381
382 /************************************************************
383 *
384 * The following routines are the exported ones
385 *
386 ************************************************************/
387
read_trr_natoms(char * fn,int * natoms)388 int read_trr_natoms(char *fn, int *natoms) {
389 XDRFILE *xd;
390 t_trnheader sh;
391 int result;
392
393 xd = xdrfile_open(fn, "r");
394 if (NULL == xd)
395 return exdrFILENOTFOUND;
396 if ((result = do_trnheader(xd, 1, &sh)) != exdrOK)
397 return result;
398 xdrfile_close(xd);
399 *natoms = sh.natoms;
400
401 return exdrOK;
402 }
403
404
write_trr(XDRFILE * xd,int natoms,int step,float t,float lambda,matrix box,rvec * x,rvec * v,rvec * f)405 int write_trr(XDRFILE *xd, int natoms, int step, float t, float lambda,
406 matrix box, rvec *x, rvec *v, rvec *f) {
407 int *plcholder;
408 return do_trn(xd, 0, &step, &t, &lambda, box, &natoms, x, v, f, plcholder);
409 }
410
read_trr(XDRFILE * xd,int natoms,int * step,float * t,float * lambda,matrix box,rvec * x,rvec * v,rvec * f,int * has_prop)411 int read_trr(XDRFILE *xd, int natoms, int *step, float *t, float *lambda,
412 matrix box, rvec *x, rvec *v, rvec *f, int *has_prop) {
413 return do_trn(xd, 1, step, t, lambda, box, &natoms, x, v, f, has_prop);
414 }
415