Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2009-2014, Erik Lindahl & David van der Spoel
3 : All rights reserved.
4 :
5 : Redistribution and use in source and binary forms, with or without
6 : modification, are permitted provided that the following conditions are met:
7 :
8 : 1. Redistributions of source code must retain the above copyright notice, this
9 : list of conditions and the following disclaimer.
10 :
11 : 2. Redistributions in binary form must reproduce the above copyright notice,
12 : this list of conditions and the following disclaimer in the documentation
13 : and/or other materials provided with the distribution.
14 :
15 : THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
16 : AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
17 : IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
18 : DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
19 : FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
20 : DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
21 : SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
22 : CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
23 : OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
24 : OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
26 : /* -*- mode: c; tab-width: 4; indent-tabs-mode: t; c-basic-offset: 4 -*-
27 : *
28 : * $Id$
29 : *
30 : /* -*- mode: c; tab-width: 4; indent-tabs-mode: t; c-basic-offset: 4 -*-
31 : *
32 : * $Id$
33 : *
34 : * Copyright (c) 2009-2014, Erik Lindahl & David van der Spoel
35 : * All rights reserved.
36 : *
37 : * Redistribution and use in source and binary forms, with or without
38 : * modification, are permitted provided that the following conditions are met:
39 : *
40 : * 1. Redistributions of source code must retain the above copyright notice, this
41 : * list of conditions and the following disclaimer.
42 : *
43 : * 2. Redistributions in binary form must reproduce the above copyright notice,
44 : * this list of conditions and the following disclaimer in the documentation
45 : * and/or other materials provided with the distribution.
46 : *
47 : * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
48 : * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
49 : * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
50 : * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
51 : * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
52 : * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
53 : * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
54 : * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
55 : * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
56 : * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
57 : */
58 :
59 : #include <stdlib.h>
60 : #include <string.h>
61 :
62 : #ifdef HAVE_CONFIG_H
63 : #include "config.h"
64 : #endif
65 :
66 : #include "xdrfile.h"
67 : #include "xdrfile_trr.h"
68 :
69 : #define BUFSIZE 128
70 : #define GROMACS_MAGIC 1993
71 :
72 : namespace PLMD {
73 : namespace xdrfile {
74 :
75 : typedef struct /* This struct describes the order and the */
76 : /* sizes of the structs in a trjfile, sizes are given in bytes. */
77 : {
78 : mybool bDouble; /* Double precision? */
79 : int ir_size; /* Backward compatibility */
80 : int e_size; /* Backward compatibility */
81 : int box_size; /* Non zero if a box is present */
82 : int vir_size; /* Backward compatibility */
83 : int pres_size; /* Backward compatibility */
84 : int top_size; /* Backward compatibility */
85 : int sym_size; /* Backward compatibility */
86 : int x_size; /* Non zero if coordinates are present */
87 : int v_size; /* Non zero if velocities are present */
88 : int f_size; /* Non zero if forces are present */
89 :
90 : int natoms; /* The total number of atoms */
91 : int step; /* Current step number */
92 : int nre; /* Backward compatibility */
93 : float tf; /* Current time */
94 : float lambdaf; /* Current value of lambda */
95 : double td; /* Current time */
96 : double lambdad; /* Current value of lambda */
97 : } t_trnheader;
98 :
99 123 : static int nFloatSize(t_trnheader *sh,int *nflsz)
100 : {
101 : int nflsize=0;
102 :
103 123 : if (sh->box_size)
104 123 : nflsize = sh->box_size/(DIM*DIM);
105 0 : else if (sh->x_size)
106 0 : nflsize = sh->x_size/(sh->natoms*DIM);
107 0 : else if (sh->v_size)
108 0 : nflsize = sh->v_size/(sh->natoms*DIM);
109 0 : else if (sh->f_size)
110 0 : nflsize = sh->f_size/(sh->natoms*DIM);
111 : else
112 : return exdrHEADER;
113 :
114 123 : if (((nflsize != sizeof(float)) && (nflsize != sizeof(double))))
115 : return exdrHEADER;
116 :
117 123 : *nflsz = nflsize;
118 :
119 123 : return exdrOK;
120 : }
121 :
122 132 : static int do_trnheader(XDRFILE *xd,mybool bRead,t_trnheader *sh)
123 : {
124 132 : int magic=GROMACS_MAGIC;
125 : int nflsz,slen,result;
126 : const char *version = "GMX_trn_file";
127 : char buf[BUFSIZE];
128 :
129 132 : if (xdrfile_read_int(&magic,1,xd) != 1)
130 : return exdrINT;
131 :
132 123 : if (bRead)
133 : {
134 75 : if (xdrfile_read_int(&slen,1,xd) != 1)
135 : return exdrINT;
136 75 : if (slen != strlen(version)+1)
137 : return exdrSTRING;
138 75 : if (xdrfile_read_string(buf,BUFSIZE,xd) <= 0)
139 : return exdrSTRING;
140 : }
141 : else
142 : {
143 48 : slen = strlen(version)+1;
144 48 : if (xdrfile_read_int(&slen,1,xd) != 1)
145 : return exdrINT;
146 48 : if (xdrfile_write_string(const_cast<char*>(version),xd) != (strlen(version)+1) )
147 : return exdrSTRING;
148 : }
149 123 : if (xdrfile_read_int(&sh->ir_size,1,xd) != 1)
150 : return exdrINT;
151 123 : if (xdrfile_read_int(&sh->e_size,1,xd) != 1)
152 : return exdrINT;
153 123 : if (xdrfile_read_int(&sh->box_size,1,xd) != 1)
154 : return exdrINT;
155 123 : if (xdrfile_read_int(&sh->vir_size,1,xd) != 1)
156 : return exdrINT;
157 123 : if (xdrfile_read_int(&sh->pres_size,1,xd) != 1)
158 : return exdrINT;
159 123 : if (xdrfile_read_int(&sh->top_size,1,xd) != 1)
160 : return exdrINT;
161 123 : if (xdrfile_read_int(&sh->sym_size,1,xd) != 1)
162 : return exdrINT;
163 123 : if (xdrfile_read_int(&sh->x_size,1,xd) != 1)
164 : return exdrINT;
165 123 : if (xdrfile_read_int(&sh->v_size,1,xd) != 1)
166 : return exdrINT;
167 123 : if (xdrfile_read_int(&sh->f_size,1,xd) != 1)
168 : return exdrINT;
169 123 : if (xdrfile_read_int(&sh->natoms,1,xd) != 1)
170 : return exdrINT;
171 :
172 123 : if ((result = nFloatSize(sh,&nflsz)) != exdrOK)
173 : return result;
174 123 : sh->bDouble = (nflsz == sizeof(double));
175 :
176 123 : if (xdrfile_read_int(&sh->step,1,xd) != 1)
177 : return exdrINT;
178 123 : if (xdrfile_read_int(&sh->nre,1,xd) != 1)
179 : return exdrINT;
180 123 : if (sh->bDouble)
181 : {
182 0 : if (xdrfile_read_double(&sh->td,1,xd) != 1)
183 : return exdrDOUBLE;
184 0 : sh->tf = sh->td;
185 0 : if (xdrfile_read_double(&sh->lambdad,1,xd) != 1)
186 : return exdrDOUBLE;
187 0 : sh->lambdaf = sh->lambdad;
188 : }
189 : else
190 : {
191 123 : if (xdrfile_read_float(&sh->tf,1,xd) != 1)
192 : return exdrFLOAT;
193 123 : sh->td = sh->tf;
194 123 : if (xdrfile_read_float(&sh->lambdaf,1,xd) != 1)
195 : return exdrFLOAT;
196 123 : sh->lambdad = sh->lambdaf;
197 : }
198 :
199 : return exdrOK;
200 : }
201 :
202 114 : static int do_htrn(XDRFILE *xd,mybool bRead,t_trnheader *sh,
203 : matrix box,rvec *x,rvec *v,rvec *f)
204 : {
205 : double pvd[DIM*DIM];
206 : double *dx=NULL;
207 : float pvf[DIM*DIM];
208 : float *fx=NULL;
209 : int i,j;
210 :
211 114 : if (sh->bDouble)
212 : {
213 0 : if (sh->box_size != 0)
214 : {
215 0 : if (!bRead)
216 : {
217 0 : for(i=0; (i<DIM); i++)
218 0 : for(j=0; (j<DIM); j++)
219 0 : if (NULL != box)
220 : {
221 0 : pvd[i*DIM+j] = box[i][j];
222 : }
223 : }
224 0 : if (xdrfile_read_double(pvd,DIM*DIM,xd) == DIM*DIM)
225 : {
226 0 : for(i=0; (i<DIM); i++)
227 0 : for(j=0; (j<DIM); j++)
228 0 : if (NULL != box)
229 : {
230 0 : box[i][j] = pvd[i*DIM+j];
231 : }
232 : }
233 : else
234 : return exdrDOUBLE;
235 : }
236 :
237 0 : if (sh->vir_size != 0)
238 : {
239 0 : if (xdrfile_read_double(pvd,DIM*DIM,xd) != DIM*DIM)
240 : return exdrDOUBLE;
241 : }
242 :
243 0 : if (sh->pres_size!= 0)
244 : {
245 0 : if (xdrfile_read_double(pvd,DIM*DIM,xd) != DIM*DIM)
246 : return exdrDOUBLE;
247 : }
248 :
249 0 : if ((sh->x_size != 0) || (sh->v_size != 0) || (sh->f_size != 0)) {
250 0 : dx = (double *)calloc(sh->natoms*DIM,sizeof(dx[0]));
251 0 : if (NULL == dx)
252 : return exdrNOMEM;
253 : }
254 0 : if (sh->x_size != 0)
255 : {
256 0 : if (!bRead)
257 : {
258 0 : for(i=0; (i<sh->natoms); i++)
259 0 : for(j=0; (j<DIM); j++)
260 0 : if (NULL != x)
261 : {
262 0 : dx[i*DIM+j] = x[i][j];
263 : }
264 : }
265 0 : if (xdrfile_read_double(dx,sh->natoms*DIM,xd) == sh->natoms*DIM)
266 : {
267 0 : if (bRead)
268 : {
269 0 : for(i=0; (i<sh->natoms); i++)
270 0 : for(j=0; (j<DIM); j++)
271 0 : if (NULL != x)
272 : {
273 0 : x[i][j] = dx[i*DIM+j];
274 : }
275 : }
276 : }
277 : else
278 : return exdrDOUBLE;
279 : }
280 0 : if (sh->v_size != 0)
281 : {
282 0 : if (!bRead)
283 : {
284 0 : for(i=0; (i<sh->natoms); i++)
285 0 : for(j=0; (j<DIM); j++)
286 0 : if (NULL != x)
287 : {
288 0 : dx[i*DIM+j] = v[i][j];
289 : }
290 : }
291 0 : if (xdrfile_read_double(dx,sh->natoms*DIM,xd) == sh->natoms*DIM)
292 : {
293 0 : for(i=0; (i<sh->natoms); i++)
294 0 : for(j=0; (j<DIM); j++)
295 0 : if (NULL != v)
296 : {
297 0 : v[i][j] = dx[i*DIM+j];
298 : }
299 : }
300 : else
301 : return exdrDOUBLE;
302 : }
303 0 : if (sh->f_size != 0)
304 : {
305 0 : if (!bRead)
306 : {
307 0 : for(i=0; (i<sh->natoms); i++)
308 0 : for(j=0; (j<DIM); j++)
309 0 : if (NULL != x)
310 : {
311 0 : dx[i*DIM+j] = f[i][j];
312 : }
313 : }
314 0 : if (xdrfile_read_double(dx,sh->natoms*DIM,xd) == sh->natoms*DIM)
315 : {
316 0 : for(i=0; (i<sh->natoms); i++)
317 : {
318 0 : for(j=0; (j<DIM); j++)
319 : {
320 0 : if (NULL != f)
321 : {
322 0 : f[i][j] = dx[i*DIM+j];
323 : }
324 : }
325 : }
326 : }
327 : else
328 : return exdrDOUBLE;
329 : }
330 0 : if ((sh->x_size != 0) || (sh->v_size != 0) || (sh->f_size != 0)) {
331 0 : free(dx);
332 : }
333 : }
334 : else
335 : /* Float */
336 : {
337 114 : if (sh->box_size != 0)
338 : {
339 114 : if (!bRead)
340 : {
341 192 : for(i=0; (i<DIM); i++)
342 576 : for(j=0; (j<DIM); j++)
343 432 : if (NULL != box)
344 : {
345 432 : pvf[i*DIM+j] = box[i][j];
346 : }
347 : }
348 114 : if (xdrfile_read_float(pvf,DIM*DIM,xd) == DIM*DIM)
349 : {
350 456 : for(i=0; (i<DIM); i++)
351 : {
352 1368 : for(j=0; (j<DIM); j++)
353 : {
354 1026 : if (NULL != box)
355 : {
356 1026 : box[i][j] = pvf[i*DIM+j];
357 : }
358 : }
359 : }
360 : }
361 : else
362 : return exdrFLOAT;
363 : }
364 :
365 114 : if (sh->vir_size != 0)
366 : {
367 0 : if (xdrfile_read_float(pvf,DIM*DIM,xd) != DIM*DIM)
368 : return exdrFLOAT;
369 : }
370 :
371 114 : if (sh->pres_size!= 0)
372 : {
373 0 : if (xdrfile_read_float(pvf,DIM*DIM,xd) != DIM*DIM)
374 : return exdrFLOAT;
375 : }
376 :
377 114 : if ((sh->x_size != 0) || (sh->v_size != 0) || (sh->f_size != 0)) {
378 114 : fx = (float *)calloc(sh->natoms*DIM,sizeof(fx[0]));
379 114 : if (NULL == fx)
380 : return exdrNOMEM;
381 : }
382 114 : if (sh->x_size != 0)
383 : {
384 114 : if (!bRead)
385 : {
386 2424 : for(i=0; (i<sh->natoms); i++)
387 9504 : for(j=0; (j<DIM); j++)
388 7128 : if (NULL != x)
389 : {
390 7128 : fx[i*DIM+j] = x[i][j];
391 : }
392 : }
393 114 : if (xdrfile_read_float(fx,sh->natoms*DIM,xd) == sh->natoms*DIM)
394 : {
395 114 : if (bRead)
396 : {
397 4728 : for(i=0; (i<sh->natoms); i++)
398 18648 : for(j=0; (j<DIM); j++)
399 13986 : if (NULL != x)
400 13986 : x[i][j] = fx[i*DIM+j];
401 : }
402 : }
403 : else
404 : return exdrFLOAT;
405 : }
406 114 : if (sh->v_size != 0)
407 : {
408 0 : if (!bRead)
409 : {
410 0 : for(i=0; (i<sh->natoms); i++)
411 0 : for(j=0; (j<DIM); j++)
412 0 : if (NULL != x)
413 : {
414 0 : fx[i*DIM+j] = v[i][j];
415 : }
416 : }
417 0 : if (xdrfile_read_float(fx,sh->natoms*DIM,xd) == sh->natoms*DIM)
418 : {
419 0 : for(i=0; (i<sh->natoms); i++)
420 0 : for(j=0; (j<DIM); j++)
421 0 : if (NULL != v)
422 0 : v[i][j] = fx[i*DIM+j];
423 : }
424 : else
425 : return exdrFLOAT;
426 : }
427 114 : if (sh->f_size != 0)
428 : {
429 0 : if (!bRead)
430 : {
431 0 : for(i=0; (i<sh->natoms); i++)
432 0 : for(j=0; (j<DIM); j++)
433 0 : if (NULL != x)
434 : {
435 0 : fx[i*DIM+j] = f[i][j];
436 : }
437 : }
438 0 : if (xdrfile_read_float(fx,sh->natoms*DIM,xd) == sh->natoms*DIM)
439 : {
440 0 : for(i=0; (i<sh->natoms); i++)
441 0 : for(j=0; (j<DIM); j++)
442 0 : if (NULL != f)
443 0 : f[i][j] = fx[i*DIM+j];
444 : }
445 : else
446 : return exdrFLOAT;
447 : }
448 114 : if ((sh->x_size != 0) || (sh->v_size != 0) || (sh->f_size != 0)) {
449 114 : free(fx);
450 : }
451 : }
452 : return exdrOK;
453 : }
454 :
455 123 : static int do_trn(XDRFILE *xd,mybool bRead,int *step,float *t,float *lambda,
456 : matrix box,int *natoms,rvec *x,rvec *v,rvec *f)
457 : {
458 : t_trnheader *sh;
459 : int result;
460 :
461 123 : sh = (t_trnheader *)calloc(1,sizeof(*sh));
462 :
463 123 : if (!bRead) {
464 48 : sh->box_size = (NULL != box) ? sizeof(matrix):0;
465 48 : sh->x_size = ((NULL != x) ? (*natoms*sizeof(x[0])):0);
466 48 : sh->v_size = ((NULL != v) ? (*natoms*sizeof(v[0])):0);
467 48 : sh->f_size = ((NULL != f) ? (*natoms*sizeof(f[0])):0);
468 48 : sh->natoms = *natoms;
469 48 : sh->step = *step;
470 48 : sh->nre = 0;
471 48 : sh->td = *t;
472 48 : sh->lambdad = *lambda;
473 48 : sh->tf = *t;
474 48 : sh->lambdaf = *lambda;
475 : }
476 123 : if ((result = do_trnheader(xd,bRead,sh)) != exdrOK)
477 : return result;
478 114 : if (bRead) {
479 66 : *natoms = sh->natoms;
480 66 : *step = sh->step;
481 66 : *t = sh->td;
482 66 : *lambda = sh->lambdad;
483 : }
484 114 : if ((result = do_htrn(xd,bRead,sh,box,x,v,f)) != exdrOK)
485 : return result;
486 :
487 114 : free(sh);
488 :
489 114 : return exdrOK;
490 : }
491 :
492 : /************************************************************
493 : *
494 : * The following routines are the exported ones
495 : *
496 : ************************************************************/
497 :
498 9 : int read_trr_natoms(char *fn,int *natoms)
499 : {
500 : XDRFILE *xd;
501 : t_trnheader sh;
502 : int result;
503 :
504 9 : xd = xdrfile_open(fn,"r");
505 9 : if (NULL == xd)
506 : return exdrFILENOTFOUND;
507 9 : if ((result = do_trnheader(xd,1,&sh)) != exdrOK)
508 : return result;
509 9 : xdrfile_close(xd);
510 9 : *natoms = sh.natoms;
511 :
512 9 : return exdrOK;
513 : }
514 :
515 48 : int write_trr(XDRFILE *xd,int natoms,int step,float t,float lambda,
516 : matrix box,rvec *x,rvec *v,rvec *f)
517 : {
518 48 : return do_trn(xd,0,&step,&t,&lambda,box,&natoms,x,v,f);
519 : }
520 :
521 75 : int read_trr(XDRFILE *xd,int natoms,int *step,float *t,float *lambda,
522 : matrix box,rvec *x,rvec *v,rvec *f)
523 : {
524 75 : return do_trn(xd,1,step,t,lambda,box,&natoms,x,v,f);
525 : }
526 :
527 : }
528 : }
529 :
|
