src/contrib/anaf.c

   1 /*
   2  * $Id$
   3  *
   4  *       This source code is part of
   5  *
   6  *        G   R   O   M   A   C   S
   7  *
   8  * GROningen MAchine for Chemical Simulations
   9  *
  10  *               VERSION 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * Great Red Oystrich Makes All Chemists Sane
  28  */
  29 static char *SRCID_anaf_c = "$Id$";
  30
  31 #include <stdio.h>
  32 #include <string.h>
  33 #include <math.h>
  34 #include "main.h"
  35 #include "macros.h"
  36 #include "futil.h"
  37 #include "statutil.h"
  38 #include "copyrite.h"
  39 #include "sysstuff.h"
  40 #include "txtdump.h"
  41 #include "fatal.h"
  42 #include "xtcio.h"
  43 #include "enxio.h"
  44 #include "assert.h"
  45 #include "smalloc.h"
  46 #include "gmxfio.h"
  47 #include "tpxio.h"
  48 #include "trnio.h"
  49 #include "txtdump.h"
  50 #include "vec.h"
  51
  52 static char *nm[5]  = { "OW", "HW1", "HW2", "DW", "SW" };
  53
  54 static void list_trn(char *fn)
  55 {
  56   static real mass[5] = { 15.9994, 1.008, 1.008, 0.0, 0.0 };
  57   int         i,j,m,fpread,fpwrite,nframe;
  58   rvec        *x,*v,*f,fmol[2],xcm[2],torque[j],dx;
  59   real        mmm,len;
  60   matrix      box;
  61   t_trnheader trn;
  62   bool        bOK;
  63
  64   printf("Going to open %s\n",fn);
  65   fpread  = open_trn(fn,"r");
  66   fpwrite = open_tpx(NULL,"w");
  67   fio_setdebug(fpwrite,TRUE);
  68
  69   mmm=mass[0]+2*mass[1];
  70   for(i=0; (i<5); i++)
  71     mass[i] /= mmm;
  72
  73   nframe = 0;
  74   while (fread_trnheader(fpread,&trn,&bOK)) {
  75     snew(x,trn.natoms);
  76     snew(v,trn.natoms);
  77     snew(f,trn.natoms);
  78     if (fread_htrn(fpread,&trn,
  79                    trn.box_size ? box : NULL,
  80                    trn.x_size   ? x : NULL,
  81                    trn.v_size   ? v : NULL,
  82                    trn.f_size   ? f : NULL)) {
  83
  84       if (trn.x_size && trn.f_size) {
  85         printf("There are %d atoms\n",trn.natoms);
  86         for(j=0; (j<2); j++) {
  87           clear_rvec(xcm[j]);
  88           clear_rvec(fmol[j]);
  89           clear_rvec(torque[j]);
  90           for(i=5*j; (i<5*j+5); i++) {
  91             rvec_inc(fmol[j],f[i]);
  92             for(m=0; (m<DIM); m++)
  93               xcm[j][m] += mass[i%5]*x[i][m];
  94           }
  95           for(i=5*j; (i<5*j+5); i++) {
  96             rvec_dec(x[i],xcm[j]);
  97             oprod(x[i],f[i],dx);
  98             rvec_inc(torque[j],dx);
  99             rvec_inc(x[i],xcm[j]);
 100           }
 101         }
 102         pr_rvecs(stdout,0,"FMOL  ",fmol,2);
 103         pr_rvecs(stdout,0,"TORQUE",torque,2);
 104         printf("Distance matrix Water1-Water2\n%5s","");
 105         for(j=0; (j<5); j++)
 106           printf("  %10s",nm[j]);
 107         printf("\n");
 108         for(j=0; (j<5); j++) {
 109           printf("%5s",nm[j]);
 110           for(i=5; (i<10); i++) {
 111             rvec_sub(x[i],x[j],dx);
 112             len = sqrt(iprod(dx,dx));
 113             printf("  %10.7f",len);
 114           }
 115           printf("\n");
 116         }
 117       }
 118     }
 119     sfree(x);
 120     sfree(v);
 121     sfree(f);
 122     nframe++;
 123   }
 124   if (!bOK)
 125     fprintf(stderr,"\nWARNING: Incomplete frame header: nr %d, t=%g\n",
 126             nframe,trn.t);
 127   close_tpx(fpwrite);
 128   close_trn(fpread);
 129 }
 130
 131 int main(int argc,char *argv[])
 132 {
 133   static char *desc[] = {
 134     "gmxdump reads a run input file ([TT].tpa[tt]/[TT].tpr[tt]/[TT].tpb[tt]),",
 135     "a trajectory ([TT].trj[tt]/[TT].trr[tt]/[TT].xtc[tt]) or an energy",
 136     "file ([TT].ene[tt]/[TT].edr[tt]) and prints that to standard",
 137     "output in a readable format. This program is essential for",
 138     "checking your run input file in case of problems.[PAR]"
 139   };
 140   t_filenm fnm[] = {
 141     { efTRN, "-f", NULL, ffOPTRD }
 142   };
 143 #define NFILE asize(fnm)
 144   char *fn;
 145
 146   /* Command line options */
 147
 148   CopyRight(stdout,argv[0]);
 149   parse_common_args(&argc,argv,0,FALSE,NFILE,fnm,0,NULL,
 150                     asize(desc),desc,0,NULL);
 151
 152   if (ftp2bSet(efTRN,NFILE,fnm)) {
 153     fn = ftp2fn(efTRN,NFILE,fnm);
 154     printf("Going to open %s\n",fn);
 155     list_trn(fn);
 156   }
 157
 158   thanx(stderr);
 159
 160   return 0;
 161 }