src/xdrfile_trr.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  * $Id$
   4  *
   5  * Copyright (c) Erik Lindahl, David van der Spoel 2003,2004.
   6  * Coordinate compression (c) by Frans van Hoesel.
   7  *
   8  * IN contrast to the rest of Gromacs, XDRFILE is distributed under the
   9  * BSD license, so you can use it any way you wish, including closed source:
  10  *
  11  * Permission is hereby granted, free of charge, to any person obtaining a
  12  * copy of this software and associated documentation files (the "Software"),
  13  * to deal in the Software without restriction, including without limitation
  14  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
  15  * and/or sell copies of the Software, and to permit persons to whom the
  16  * Software is furnished to do so, subject to the following conditions:
  17  *
  18  * The above copyright notice and this permission notice shall be included in
  19  * all copies or substantial portions of the Software.
  20  *
  21  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  22  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  23  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  24  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  25  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  26  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  27  * DEALINGS IN THE SOFTWARE.
  28  */
  29
  30 #include <stdlib.h>
  31 #include <string.h>
  32
  33 #ifdef HAVE_CONFIG_H
  34 #include "config.h"
  35 #endif
  36
  37 #include "xdrfile.h"
  38 #include "xdrfile_trr.h"
  39
  40 #define BUFSIZE         128
  41 #define GROMACS_MAGIC   1993
  42
  43 typedef struct          /* This struct describes the order and the      */
  44 /* sizes of the structs in a trjfile, sizes are given in bytes. */
  45 {
  46     bool  bDouble;        /* Double precision?                            */
  47     int ir_size;        /* Backward compatibility                       */
  48     int e_size;         /* Backward compatibility                       */
  49     int box_size;       /* Non zero if a box is present                 */
  50     int   vir_size;       /* Backward compatibility                     */
  51     int   pres_size;      /* Backward compatibility                     */
  52     int top_size;       /* Backward compatibility                       */
  53     int sym_size;       /* Backward compatibility                       */
  54     int x_size;         /* Non zero if coordinates are present          */
  55     int v_size;         /* Non zero if velocities are present           */
  56     int f_size;         /* Non zero if forces are present               */
  57
  58     int natoms;         /* The total number of atoms                    */
  59     int step;           /* Current step number                          */
  60     int nre;            /* Backward compatibility                       */
  61     float       tf;             /* Current time                                 */
  62     float       lambdaf;                /* Current value of lambda                      */
  63     double      td;             /* Current time                                 */
  64     double      lambdad;                /* Current value of lambda                      */
  65 } t_trnheader;
  66
  67 static int nFloatSize(t_trnheader *sh,int *nflsz)
  68 {
  69     int nflsize=0;
  70
  71     if (sh->box_size)
  72         nflsize = sh->box_size/(DIM*DIM);
  73     else if (sh->x_size)
  74         nflsize = sh->x_size/(sh->natoms*DIM);
  75     else if (sh->v_size)
  76         nflsize = sh->v_size/(sh->natoms*DIM);
  77     else if (sh->f_size)
  78         nflsize = sh->f_size/(sh->natoms*DIM);
  79     else
  80         return exdrHEADER;
  81
  82     if (((nflsize != sizeof(float)) && (nflsize != sizeof(double))))
  83         return exdrHEADER;
  84
  85     *nflsz = nflsize;
  86
  87     return exdrOK;
  88 }
  89
  90 static int do_trnheader(XDRFILE *xd,bool bRead,t_trnheader *sh)
  91 {
  92         int magic=GROMACS_MAGIC;
  93         int nflsz,slen,result;
  94         char *version = "GMX_trn_file";
  95         char buf[BUFSIZE];
  96
  97         if (xdrfile_read_int(&magic,1,xd) != 1)
  98                 return exdrINT;
  99
 100         if (bRead)
 101     {
 102         if (xdrfile_read_int(&slen,1,xd) != 1)
 103             return exdrINT;
 104         if (slen != strlen(version)+1)
 105             return exdrSTRING;
 106         if (xdrfile_read_string(buf,BUFSIZE,xd) <= 0)
 107             return exdrSTRING;
 108     }
 109         else
 110     {
 111         slen = strlen(version)+1;
 112         if (xdrfile_read_int(&slen,1,xd) != 1)
 113             return exdrINT;
 114         if (xdrfile_write_string(version,xd) != (strlen(version)+1) )
 115             return exdrSTRING;
 116     }
 117         if (xdrfile_read_int(&sh->ir_size,1,xd) != 1)
 118                 return exdrINT;
 119         if (xdrfile_read_int(&sh->e_size,1,xd) != 1)
 120                 return exdrINT;
 121         if (xdrfile_read_int(&sh->box_size,1,xd) != 1)
 122                 return exdrINT;
 123         if (xdrfile_read_int(&sh->vir_size,1,xd) != 1)
 124                 return exdrINT;
 125         if (xdrfile_read_int(&sh->pres_size,1,xd) != 1)
 126                 return exdrINT;
 127         if (xdrfile_read_int(&sh->top_size,1,xd) != 1)
 128                 return exdrINT;
 129         if (xdrfile_read_int(&sh->sym_size,1,xd) != 1)
 130                 return exdrINT;
 131         if (xdrfile_read_int(&sh->x_size,1,xd) != 1)
 132                 return exdrINT;
 133         if (xdrfile_read_int(&sh->v_size,1,xd) != 1)
 134                 return exdrINT;
 135         if (xdrfile_read_int(&sh->f_size,1,xd) != 1)
 136                 return exdrINT;
 137         if (xdrfile_read_int(&sh->natoms,1,xd) != 1)
 138                 return exdrINT;
 139
 140         if ((result = nFloatSize(sh,&nflsz)) != exdrOK)
 141                 return result;
 142         sh->bDouble = (nflsz == sizeof(double));
 143
 144         if (xdrfile_read_int(&sh->step,1,xd) != 1)
 145                 return exdrINT;
 146         if (xdrfile_read_int(&sh->nre,1,xd) != 1)
 147                 return exdrINT;
 148         if (sh->bDouble)
 149     {
 150         if (xdrfile_read_double(&sh->td,1,xd) != 1)
 151             return exdrDOUBLE;
 152         sh->tf = sh->td;
 153         if (xdrfile_read_double(&sh->lambdad,1,xd) != 1)
 154             return exdrDOUBLE;
 155         sh->lambdaf = sh->lambdad;
 156     }
 157         else
 158     {
 159         if (xdrfile_read_float(&sh->tf,1,xd) != 1)
 160             return exdrFLOAT;
 161         sh->td = sh->tf;
 162         if (xdrfile_read_float(&sh->lambdaf,1,xd) != 1)
 163             return exdrFLOAT;
 164         sh->lambdad = sh->lambdaf;
 165     }
 166
 167     return exdrOK;
 168 }
 169
 170 static int do_htrn(XDRFILE *xd,bool bRead,t_trnheader *sh,
 171                                    matrix box,rvec *x,rvec *v,rvec *f)
 172 {
 173         double pvd[DIM*DIM];
 174         double *dx;
 175         float  pvf[DIM*DIM];
 176         float  *fx;
 177         int    i,j;
 178
 179         if (sh->bDouble)
 180         {
 181                 if (sh->box_size != 0)
 182         {
 183             if (!bRead)
 184             {
 185                 for(i=0; (i<DIM); i++)
 186                     for(j=0; (j<DIM); j++)
 187                         if (NULL != box)
 188                         {
 189                             pvd[i*DIM+j] = box[i][j];
 190                         }
 191             }
 192             if (xdrfile_read_double(pvd,DIM*DIM,xd) == DIM*DIM)
 193             {
 194                 for(i=0; (i<DIM); i++)
 195                     for(j=0; (j<DIM); j++)
 196                         if (NULL != box)
 197                         {
 198                             box[i][j] = pvd[i*DIM+j];
 199                         }
 200             }
 201             else
 202                 return exdrDOUBLE;
 203         }
 204
 205                 if (sh->vir_size != 0)
 206         {
 207             if (xdrfile_read_double(pvd,DIM*DIM,xd) != DIM*DIM)
 208                 return exdrDOUBLE;
 209         }
 210
 211                 if (sh->pres_size!= 0)
 212         {
 213             if (xdrfile_read_double(pvd,DIM*DIM,xd) != DIM*DIM)
 214                 return exdrDOUBLE;
 215         }
 216
 217                 if ((sh->x_size != 0) || (sh->v_size != 0) || (sh->f_size != 0)) {
 218                         dx = calloc(sh->natoms*DIM,sizeof(dx[0]));
 219                         if (NULL == dx)
 220                                 return exdrNOMEM;
 221                 }
 222                 if (sh->x_size   != 0)
 223         {
 224             if (!bRead)
 225             {
 226                 for(i=0; (i<sh->natoms); i++)
 227                     for(j=0; (j<DIM); j++)
 228                         if (NULL != x)
 229                         {
 230                             dx[i*DIM+j] = x[i][j];
 231                         }
 232             }
 233             if (xdrfile_read_double(dx,sh->natoms*DIM,xd) == sh->natoms*DIM)
 234             {
 235                 if (bRead)
 236                 {
 237                     for(i=0; (i<sh->natoms); i++)
 238                         for(j=0; (j<DIM); j++)
 239                             if (NULL != x)
 240                             {
 241                                 x[i][j] = dx[i*DIM+j];
 242                             }
 243                 }
 244             }
 245             else
 246                 return exdrDOUBLE;
 247         }
 248                 if (sh->v_size   != 0)
 249         {
 250             if (!bRead)
 251             {
 252                 for(i=0; (i<sh->natoms); i++)
 253                     for(j=0; (j<DIM); j++)
 254                         if (NULL != x)
 255                         {
 256                             dx[i*DIM+j] = v[i][j];
 257                         }
 258             }
 259             if (xdrfile_read_double(dx,sh->natoms*DIM,xd) == sh->natoms*DIM)
 260             {
 261                 for(i=0; (i<sh->natoms); i++)
 262                     for(j=0; (j<DIM); j++)
 263                         if (NULL != v)
 264                         {
 265                             v[i][j] = dx[i*DIM+j];
 266                         }
 267             }
 268             else
 269                 return exdrDOUBLE;
 270         }
 271                 if (sh->f_size   != 0)
 272         {
 273             if (!bRead)
 274             {
 275                 for(i=0; (i<sh->natoms); i++)
 276                     for(j=0; (j<DIM); j++)
 277                         if (NULL != x)
 278                         {
 279                             dx[i*DIM+j] = f[i][j];
 280                         }
 281             }
 282             if (xdrfile_read_double(dx,sh->natoms*DIM,xd) == sh->natoms*DIM)
 283             {
 284                 for(i=0; (i<sh->natoms); i++)
 285                 {
 286                     for(j=0; (j<DIM); j++)
 287                     {
 288                         if (NULL != f)
 289                         {
 290                             f[i][j] = dx[i*DIM+j];
 291                         }
 292                     }
 293                 }
 294             }
 295             else
 296                 return exdrDOUBLE;
 297         }
 298                 if ((sh->x_size != 0) || (sh->v_size != 0) || (sh->f_size != 0)) {
 299                         free(dx);
 300                 }
 301         }
 302         else
 303                 /* Float */
 304         {
 305                 if (sh->box_size != 0)
 306         {
 307             if (!bRead)
 308             {
 309                 for(i=0; (i<DIM); i++)
 310                     for(j=0; (j<DIM); j++)
 311                         if (NULL != box)
 312                         {
 313                             pvf[i*DIM+j] = box[i][j];
 314                         }
 315             }
 316             if (xdrfile_read_float(pvf,DIM*DIM,xd) == DIM*DIM)
 317             {
 318                 for(i=0; (i<DIM); i++)
 319                 {
 320                     for(j=0; (j<DIM); j++)
 321                     {
 322                         if (NULL != box)
 323                         {
 324                             box[i][j] = pvf[i*DIM+j];
 325                         }
 326                     }
 327                 }
 328             }
 329             else
 330                 return exdrFLOAT;
 331         }
 332
 333                 if (sh->vir_size != 0)
 334         {
 335             if (xdrfile_read_float(pvf,DIM*DIM,xd) != DIM*DIM)
 336                 return exdrFLOAT;
 337         }
 338
 339                 if (sh->pres_size!= 0)
 340         {
 341             if (xdrfile_read_float(pvf,DIM*DIM,xd) != DIM*DIM)
 342                 return exdrFLOAT;
 343         }
 344
 345                 if ((sh->x_size != 0) || (sh->v_size != 0) || (sh->f_size != 0)) {
 346                         fx = calloc(sh->natoms*DIM,sizeof(fx[0]));
 347                         if (NULL == fx)
 348                                 return exdrNOMEM;
 349                 }
 350                 if (sh->x_size   != 0)
 351         {
 352             if (!bRead)
 353             {
 354                 for(i=0; (i<sh->natoms); i++)
 355                     for(j=0; (j<DIM); j++)
 356                         if (NULL != x)
 357                         {
 358                             fx[i*DIM+j] = x[i][j];
 359                         }
 360             }
 361             if (xdrfile_read_float(fx,sh->natoms*DIM,xd) == sh->natoms*DIM)
 362             {
 363                 if (bRead)
 364                 {
 365                     for(i=0; (i<sh->natoms); i++)
 366                         for(j=0; (j<DIM); j++)
 367                             if (NULL != x)
 368                                 x[i][j] = fx[i*DIM+j];
 369                 }
 370             }
 371             else
 372                 return exdrFLOAT;
 373         }
 374                 if (sh->v_size   != 0)
 375         {
 376             if (!bRead)
 377             {
 378                 for(i=0; (i<sh->natoms); i++)
 379                     for(j=0; (j<DIM); j++)
 380                         if (NULL != x)
 381                         {
 382                             fx[i*DIM+j] = v[i][j];
 383                         }
 384             }
 385             if (xdrfile_read_float(fx,sh->natoms*DIM,xd) == sh->natoms*DIM)
 386             {
 387                 for(i=0; (i<sh->natoms); i++)
 388                     for(j=0; (j<DIM); j++)
 389                         if (NULL != v)
 390                             v[i][j] = fx[i*DIM+j];
 391             }
 392             else
 393                 return exdrFLOAT;
 394         }
 395                 if (sh->f_size   != 0)
 396         {
 397            if (!bRead)
 398             {
 399                 for(i=0; (i<sh->natoms); i++)
 400                     for(j=0; (j<DIM); j++)
 401                         if (NULL != x)
 402                         {
 403                             fx[i*DIM+j] = f[i][j];
 404                         }
 405             }
 406              if (xdrfile_read_float(fx,sh->natoms*DIM,xd) == sh->natoms*DIM)
 407             {
 408                 for(i=0; (i<sh->natoms); i++)
 409                     for(j=0; (j<DIM); j++)
 410                         if (NULL != f)
 411                             f[i][j] = fx[i*DIM+j];
 412             }
 413             else
 414                 return exdrFLOAT;
 415         }
 416                 if ((sh->x_size != 0) || (sh->v_size != 0) || (sh->f_size != 0)) {
 417                         free(fx);
 418                 }
 419         }
 420         return exdrOK;
 421 }
 422
 423 static int do_trn(XDRFILE *xd,bool bRead,int *step,float *t,float *lambda,
 424                                   matrix box,int *natoms,rvec *x,rvec *v,rvec *f)
 425 {
 426     t_trnheader *sh;
 427     int result;
 428
 429     sh = calloc(1,sizeof(*sh));
 430
 431     if (!bRead) {
 432         sh->box_size = (NULL != box) ? sizeof(matrix):0;
 433         sh->x_size   = ((NULL != x) ? (*natoms*sizeof(x[0])):0);
 434         sh->v_size   = ((NULL != v) ? (*natoms*sizeof(v[0])):0);
 435         sh->f_size   = ((NULL != f) ? (*natoms*sizeof(f[0])):0);
 436         sh->natoms = *natoms;
 437         sh->step   = *step;
 438         sh->nre    = 0;
 439         sh->td      = *t;
 440         sh->lambdad = *lambda;
 441         sh->tf      = *t;
 442         sh->lambdaf = *lambda;
 443     }
 444     if ((result = do_trnheader(xd,bRead,sh)) != exdrOK)
 445         return result;
 446     if (bRead) {
 447         *natoms = sh->natoms;
 448         *step   = sh->step;
 449         *t      = sh->td;
 450         *lambda = sh->lambdad;
 451     }
 452     if ((result = do_htrn(xd,bRead,sh,box,x,v,f)) != exdrOK)
 453         return result;
 454
 455     free(sh);
 456
 457     return exdrOK;
 458 }
 459
 460 /************************************************************
 461  *
 462  *  The following routines are the exported ones
 463  *
 464  ************************************************************/
 465
 466 int read_trr_natoms(char *fn,int *natoms)
 467 {
 468         XDRFILE *xd;
 469         t_trnheader sh;
 470         int  result;
 471
 472         xd = xdrfile_open(fn,"r");
 473         if (NULL == xd)
 474                 return exdrFILENOTFOUND;
 475         if ((result = do_trnheader(xd,1,&sh)) != exdrOK)
 476                 return result;
 477         xdrfile_close(xd);
 478         *natoms = sh.natoms;
 479
 480         return exdrOK;
 481 }
 482
 483 int write_trr(XDRFILE *xd,int natoms,int step,float t,float lambda,
 484                           matrix box,rvec *x,rvec *v,rvec *f)
 485 {
 486         return do_trn(xd,0,&step,&t,&lambda,box,&natoms,x,v,f);
 487 }
 488
 489 int read_trr(XDRFILE *xd,int natoms,int *step,float *t,float *lambda,
 490                          matrix box,rvec *x,rvec *v,rvec *f)
 491 {
 492         return do_trn(xd,1,step,t,lambda,box,&natoms,x,v,f);
 493 }
 494