src/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       f = *r;
  60     ret = xdr_float(xdrs,&f);
  61     if (bRead)
  62       *r = f;
  63
  64     return ret;
  65 #else
  66     return xdr_float(xdrs,(float *)r);
  67 #endif
  68 }
  69
  70
  71 t_fileio *open_xtc(const char *fn,const char *mode)
  72 {
  73     return gmx_fio_open(fn,mode);
  74 }
  75
  76 void close_xtc(t_fileio *fio)
  77 {
  78     gmx_fio_close(fio);
  79 }
  80
  81 static void check_xtc_magic(int magic)
  82 {
  83   if (magic != XTC_MAGIC)
  84     gmx_fatal(FARGS,"Magic Number Error in XTC file (read %d, should be %d)",
  85                 magic,XTC_MAGIC);
  86 }
  87
  88 int xtc_check(const char *str,gmx_bool bResult,const char *file,int line)
  89 {
  90   if (!bResult) {
  91     if (debug)
  92       fprintf(debug,"\nXTC error: read/write of %s failed, "
  93               "source file %s, line %d\n",str,file,line);
  94     return 0;
  95   }
  96   return 1;
  97 }
  98
  99 void xtc_check_fat_err(const char *str,gmx_bool bResult,const char *file,int line)
 100 {
 101   if (!bResult) {
 102     gmx_fatal(FARGS,"XTC read/write of %s failed, "
 103                 "source file %s, line %d\n",str,file,line);
 104   }
 105 }
 106
 107 static int xtc_header(XDR *xd,int *magic,int *natoms,int *step,real *time,
 108                       gmx_bool bRead,gmx_bool *bOK)
 109 {
 110   int result;
 111
 112   if (xdr_int(xd,magic) == 0)
 113     return 0;
 114   result=XTC_CHECK("natoms", xdr_int(xd,natoms));  /* number of atoms */
 115   if (result)
 116     result=XTC_CHECK("step",   xdr_int(xd,step));    /* frame number    */
 117   if (result)
 118     result=XTC_CHECK("time",   xdr_r2f(xd,time,bRead));   /* time */
 119   *bOK=(result!=0);
 120
 121   return result;
 122 }
 123
 124 static int xtc_coord(XDR *xd,int *natoms,matrix box,rvec *x,real *prec, gmx_bool bRead)
 125 {
 126   int i,j,result;
 127 #ifdef GMX_DOUBLE
 128   float *ftmp;
 129   float fprec;
 130 #endif
 131
 132   /* box */
 133   result=1;
 134   for(i=0; ((i<DIM) && result); i++)
 135     for(j=0; ((j<DIM) && result); j++)
 136       result=XTC_CHECK("box",xdr_r2f(xd,&(box[i][j]),bRead));
 137
 138   if (!result)
 139       return result;
 140
 141 #ifdef GMX_DOUBLE
 142   /* allocate temp. single-precision array */
 143   snew(ftmp,(*natoms)*DIM);
 144
 145   /* Copy data to temp. array if writing */
 146   if(!bRead)
 147   {
 148       for(i=0; (i<*natoms); i++)
 149       {
 150           ftmp[DIM*i+XX]=x[i][XX];
 151           ftmp[DIM*i+YY]=x[i][YY];
 152           ftmp[DIM*i+ZZ]=x[i][ZZ];
 153       }
 154       fprec = *prec;
 155   }
 156   result=XTC_CHECK("x",xdr3dfcoord(xd,ftmp,natoms,&fprec));
 157
 158   /* Copy from temp. array if reading */
 159   if(bRead)
 160   {
 161       for(i=0; (i<*natoms); i++)
 162       {
 163           x[i][XX] = ftmp[DIM*i+XX];
 164           x[i][YY] = ftmp[DIM*i+YY];
 165           x[i][ZZ] = ftmp[DIM*i+ZZ];
 166       }
 167       *prec = fprec;
 168   }
 169   sfree(ftmp);
 170 #else
 171     result=XTC_CHECK("x",xdr3dfcoord(xd,x[0],natoms,prec));
 172 #endif
 173
 174   return result;
 175 }
 176
 177
 178
 179 int write_xtc(t_fileio *fio,
 180               int natoms,int step,real time,
 181               matrix box,rvec *x,real prec)
 182 {
 183   int magic_number = XTC_MAGIC;
 184   XDR *xd;
 185   gmx_bool bDum;
 186   int bOK;
 187
 188   xd = gmx_fio_getxdr(fio);
 189   /* write magic number and xtc identidier */
 190   if (xtc_header(xd,&magic_number,&natoms,&step,&time,FALSE,&bDum) == 0)
 191   {
 192           return 0;
 193   }
 194
 195   /* write data */
 196   bOK = xtc_coord(xd,&natoms,box,x,&prec,FALSE); /* bOK will be 1 if writing went well */
 197
 198   if(bOK)
 199   {
 200           if(gmx_fio_flush(fio) !=0)
 201           {
 202                   bOK = 0;
 203           }
 204   }
 205   return bOK;  /* 0 if bad, 1 if writing went well */
 206 }
 207
 208 int read_first_xtc(t_fileio *fio,int *natoms,int *step,real *time,
 209                    matrix box,rvec **x,real *prec,gmx_bool *bOK)
 210 {
 211   int magic;
 212   XDR *xd;
 213
 214   *bOK=TRUE;
 215   xd = gmx_fio_getxdr(fio);
 216
 217   /* read header and malloc x */
 218   if ( !xtc_header(xd,&magic,natoms,step,time,TRUE,bOK))
 219     return 0;
 220
 221   /* Check magic number */
 222   check_xtc_magic(magic);
 223
 224   snew(*x,*natoms);
 225
 226   *bOK=xtc_coord(xd,natoms,box,*x,prec,TRUE);
 227
 228   return *bOK;
 229 }
 230
 231 int read_next_xtc(t_fileio* fio,
 232                   int natoms,int *step,real *time,
 233                   matrix box,rvec *x,real *prec,gmx_bool *bOK)
 234 {
 235   int magic;
 236   int n;
 237   XDR *xd;
 238
 239   *bOK=TRUE;
 240   xd = gmx_fio_getxdr(fio);
 241
 242   /* read header */
 243   if (!xtc_header(xd,&magic,&n,step,time,TRUE,bOK))
 244     return 0;
 245
 246   /* Check magic number */
 247   check_xtc_magic(magic);
 248
 249   if (n > natoms) {
 250     gmx_fatal(FARGS, "Frame contains more atoms (%d) than expected (%d)",
 251               n, natoms);
 252   }
 253
 254   *bOK=xtc_coord(xd,&natoms,box,x,prec,TRUE);
 255
 256   return *bOK;
 257 }
 258
 259