src/gromacs/gmxlib/xtcio.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * GROningen Mixture of Alchemy and Childrens' Stories
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <string.h>
  40 #include "typedefs.h"
  41 #include "xdrf.h"
  42 #include "gmxfio.h"
  43 #include "xtcio.h"
  44 #include "smalloc.h"
  45 #include "vec.h"
  46 #include "futil.h"
  47 #include "gmx_fatal.h"
  48
  49 #define XTC_MAGIC 1995
  50
  51
  52 static int xdr_r2f(XDR *xdrs, real *r, gmx_bool bRead)
  53 {
  54 #ifdef GMX_DOUBLE
  55     float f;
  56     int   ret;
  57
  58     if (!bRead)
  59     {
  60         f = *r;
  61     }
  62     ret = xdr_float(xdrs, &f);
  63     if (bRead)
  64     {
  65         *r = f;
  66     }
  67
  68     return ret;
  69 #else
  70     return xdr_float(xdrs, (float *)r);
  71 #endif
  72 }
  73
  74
  75 t_fileio *open_xtc(const char *fn, const char *mode)
  76 {
  77     return gmx_fio_open(fn, mode);
  78 }
  79
  80 void close_xtc(t_fileio *fio)
  81 {
  82     gmx_fio_close(fio);
  83 }
  84
  85 static void check_xtc_magic(int magic)
  86 {
  87     if (magic != XTC_MAGIC)
  88     {
  89         gmx_fatal(FARGS, "Magic Number Error in XTC file (read %d, should be %d)",
  90                   magic, XTC_MAGIC);
  91     }
  92 }
  93
  94 int xtc_check(const char *str, gmx_bool bResult, const char *file, int line)
  95 {
  96     if (!bResult)
  97     {
  98         if (debug)
  99         {
 100             fprintf(debug, "\nXTC error: read/write of %s failed, "
 101                     "source file %s, line %d\n", str, file, line);
 102         }
 103         return 0;
 104     }
 105     return 1;
 106 }
 107
 108 void xtc_check_fat_err(const char *str, gmx_bool bResult, const char *file, int line)
 109 {
 110     if (!bResult)
 111     {
 112         gmx_fatal(FARGS, "XTC read/write of %s failed, "
 113                   "source file %s, line %d\n", str, file, line);
 114     }
 115 }
 116
 117 static int xtc_header(XDR *xd, int *magic, int *natoms, int *step, real *time,
 118                       gmx_bool bRead, gmx_bool *bOK)
 119 {
 120     int result;
 121
 122     if (xdr_int(xd, magic) == 0)
 123     {
 124         return 0;
 125     }
 126     result = XTC_CHECK("natoms", xdr_int(xd, natoms)); /* number of atoms */
 127     if (result)
 128     {
 129         result = XTC_CHECK("step",   xdr_int(xd, step)); /* frame number    */
 130     }
 131     if (result)
 132     {
 133         result = XTC_CHECK("time",   xdr_r2f(xd, time, bRead)); /* time */
 134     }
 135     *bOK = (result != 0);
 136
 137     return result;
 138 }
 139
 140 static int xtc_coord(XDR *xd, int *natoms, matrix box, rvec *x, real *prec, gmx_bool bRead)
 141 {
 142     int    i, j, result;
 143 #ifdef GMX_DOUBLE
 144     float *ftmp;
 145     float  fprec;
 146 #endif
 147
 148     /* box */
 149     result = 1;
 150     for (i = 0; ((i < DIM) && result); i++)
 151     {
 152         for (j = 0; ((j < DIM) && result); j++)
 153         {
 154             result = XTC_CHECK("box", xdr_r2f(xd, &(box[i][j]), bRead));
 155         }
 156     }
 157
 158     if (!result)
 159     {
 160         return result;
 161     }
 162
 163 #ifdef GMX_DOUBLE
 164     /* allocate temp. single-precision array */
 165     snew(ftmp, (*natoms)*DIM);
 166
 167     /* Copy data to temp. array if writing */
 168     if (!bRead)
 169     {
 170         for (i = 0; (i < *natoms); i++)
 171         {
 172             ftmp[DIM*i+XX] = x[i][XX];
 173             ftmp[DIM*i+YY] = x[i][YY];
 174             ftmp[DIM*i+ZZ] = x[i][ZZ];
 175         }
 176         fprec = *prec;
 177     }
 178     result = XTC_CHECK("x", xdr3dfcoord(xd, ftmp, natoms, &fprec));
 179
 180     /* Copy from temp. array if reading */
 181     if (bRead)
 182     {
 183         for (i = 0; (i < *natoms); i++)
 184         {
 185             x[i][XX] = ftmp[DIM*i+XX];
 186             x[i][YY] = ftmp[DIM*i+YY];
 187             x[i][ZZ] = ftmp[DIM*i+ZZ];
 188         }
 189         *prec = fprec;
 190     }
 191     sfree(ftmp);
 192 #else
 193     result = XTC_CHECK("x", xdr3dfcoord(xd, x[0], natoms, prec));
 194 #endif
 195
 196     return result;
 197 }
 198
 199
 200
 201 int write_xtc(t_fileio *fio,
 202               int natoms, int step, real time,
 203               matrix box, rvec *x, real prec)
 204 {
 205     int      magic_number = XTC_MAGIC;
 206     XDR     *xd;
 207     gmx_bool bDum;
 208     int      bOK;
 209
 210     xd = gmx_fio_getxdr(fio);
 211     /* write magic number and xtc identidier */
 212     if (xtc_header(xd, &magic_number, &natoms, &step, &time, FALSE, &bDum) == 0)
 213     {
 214         return 0;
 215     }
 216
 217     /* write data */
 218     bOK = xtc_coord(xd, &natoms, box, x, &prec, FALSE); /* bOK will be 1 if writing went well */
 219
 220     if (bOK)
 221     {
 222         if (gmx_fio_flush(fio) != 0)
 223         {
 224             bOK = 0;
 225         }
 226     }
 227     return bOK; /* 0 if bad, 1 if writing went well */
 228 }
 229
 230 int read_first_xtc(t_fileio *fio, int *natoms, int *step, real *time,
 231                    matrix box, rvec **x, real *prec, gmx_bool *bOK)
 232 {
 233     int  magic;
 234     XDR *xd;
 235
 236     *bOK = TRUE;
 237     xd   = gmx_fio_getxdr(fio);
 238
 239     /* read header and malloc x */
 240     if (!xtc_header(xd, &magic, natoms, step, time, TRUE, bOK))
 241     {
 242         return 0;
 243     }
 244
 245     /* Check magic number */
 246     check_xtc_magic(magic);
 247
 248     snew(*x, *natoms);
 249
 250     *bOK = xtc_coord(xd, natoms, box, *x, prec, TRUE);
 251
 252     return *bOK;
 253 }
 254
 255 int read_next_xtc(t_fileio* fio,
 256                   int natoms, int *step, real *time,
 257                   matrix box, rvec *x, real *prec, gmx_bool *bOK)
 258 {
 259     int  magic;
 260     int  n;
 261     XDR *xd;
 262
 263     *bOK = TRUE;
 264     xd   = gmx_fio_getxdr(fio);
 265
 266     /* read header */
 267     if (!xtc_header(xd, &magic, &n, step, time, TRUE, bOK))
 268     {
 269         return 0;
 270     }
 271
 272     /* Check magic number */
 273     check_xtc_magic(magic);
 274
 275     if (n > natoms)
 276     {
 277         gmx_fatal(FARGS, "Frame contains more atoms (%d) than expected (%d)",
 278                   n, natoms);
 279     }
 280
 281     *bOK = xtc_coord(xd, &natoms, box, x, prec, TRUE);
 282
 283     return *bOK;
 284 }