src/local/do_shift.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_do_shift_c = "$Id$";
  30
  31 #include <stdlib.h>
  32 #include "errno.h"
  33 #include "sysstuff.h"
  34 #include "typedefs.h"
  35 #include "string2.h"
  36 #include "strdb.h"
  37 #include "macros.h"
  38 #include "smalloc.h"
  39 #include "mshift.h"
  40 #include "statutil.h"
  41 #include "copyrite.h"
  42 #include "confio.h"
  43 #include "fatal.h"
  44 #include "xvgr.h"
  45 #include "gstat.h"
  46 #include "rdgroup.h"
  47 #include "pdbio.h"
  48
  49 void cat(FILE *out,char *fn,real t)
  50 {
  51   FILE *in;
  52   char *ptr,buf[256];
  53   int    anr,rnr;
  54   char   anm[24],rnm[24];
  55   double f1,f2,f3,f4,f5,f6;
  56
  57   in=ffopen(fn,"r");
  58   while ((ptr=fgets2(buf,255,in)) != NULL) {
  59     sscanf(buf,"%d%d%s%s%lf%lf%lf%lf%lf%lf",
  60            &anr,&rnr,rnm,anm,&f1,&f2,&f3,&f4,&f5,&f6);
  61     fprintf(out,"%8g  %10g  %10g  %10g  %10g  %10g  %10g  %s%d-%s%d\n",
  62             t,f6,f1,f2,f3,f4,f5,rnm,rnr,anm,anr);
  63   }
  64   /*if ((int)strlen(buf) > 0)
  65     fprintf(out,"%s\n",buf);*/
  66   fflush(out);
  67   fclose(in);
  68 }
  69
  70 int main(int argc,char *argv[])
  71 {
  72   static char *desc[] = {
  73     "do_shift reads a trajectory file and computes the chemical",
  74     "shift for each time frame (or every [BB]dt[bb] ps) by",
  75     "calling the 'total' program. If you do not have the total program,",
  76     "get it. do_shift assumes that the total executable is in",
  77     "/home/mdgroup/total/total. If that is not the case, then you should",
  78     "set an environment variable [BB]TOTAL[bb] as in: [PAR]",
  79     "[TT]setenv TOTAL /usr/local/bin/total[tt][PAR]",
  80     "where the right hand side should point to the total executable.[PAR]",
  81     "Output is printed in files shift.out where t is the time of the frame.[PAR]",
  82     "The program also needs an input file called [BB]random.dat[bb] which",
  83     "contains the random coil chemical shifts of all protons."
  84   };
  85   static real dt=0.0;
  86   t_pargs pa[] = {
  87     { "-dt", FALSE, etREAL, &dt, "Time interval between frames." }
  88   };
  89   static char *bugs[] = {
  90     "The program is very slow"
  91   };
  92   static     char *OXYGEN="O";
  93   FILE       *out,*tot,*fp;
  94   t_topology *top;
  95   t_atoms    *atoms;
  96   int        status,nres;
  97   real       t,nt;
  98   int        i,natoms,nframe=0;
  99   matrix     box;
 100   int        gnx;
 101   char       *grpnm,*randf;
 102   atom_id    *index;
 103   rvec       *x,*x_s;
 104   char       pdbfile[L_tmpnam],tmpfile[L_tmpnam];
 105   char       total[256],*dptr;
 106   t_filenm   fnm[] = {
 107     { efTRX, "-f",   NULL,     ffREAD },
 108     { efTPX, NULL,   NULL,     ffREAD },
 109     { efNDX, NULL,   NULL,     ffREAD },
 110     { efOUT, "-o",   "shift",  ffWRITE },
 111     { efDAT, "-d",   "random", ffREAD }
 112   };
 113   char *leg[] = { "shift","ring","anisCO","anisCN","sigmaE","sum" };
 114 #define NFILE asize(fnm)
 115
 116   CopyRight(stdout,argv[0]);
 117   parse_common_args(&argc,argv,PCA_CAN_TIME,TRUE,NFILE,fnm,
 118                     asize(pa),pa,asize(desc),desc,asize(bugs),bugs);
 119
 120   top=read_top(ftp2fn(efTPX,NFILE,fnm));
 121   atoms=&(top->atoms);
 122   nres=atoms->nres;
 123   for(i=0; (i<atoms->nr); i++)
 124     if ((strcmp(*atoms->atomname[i],"O1") == 0) ||
 125         (strcmp(*atoms->atomname[i],"O2") == 0) ||
 126         (strcmp(*atoms->atomname[i],"OXT") == 0) ||
 127         (strcmp(*atoms->atomname[i],"OT") == 0))
 128       atoms->atomname[i]=&OXYGEN;
 129   rd_index(ftp2fn(efNDX,NFILE,fnm),1,&gnx,&index,&grpnm);
 130
 131   snew(x_s,atoms->nr);
 132
 133   (void) tmpnam(pdbfile);
 134   (void) tmpnam(tmpfile);
 135   fprintf(stderr,"pdbfile = %s\ntmpfile = %s\n",pdbfile,tmpfile);
 136
 137   if ((dptr=getenv("TOTAL")) == NULL)
 138     dptr="/home/mdgroup/total/total";
 139   sprintf(total,"%s > /dev/null",dptr);
 140   fprintf(stderr,"total cmd='%s'\n",total);
 141   randf=ftp2fn(efDAT,NFILE,fnm);
 142
 143   natoms=read_first_x(&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
 144   if (natoms != atoms->nr)
 145     fatal_error(0,"Trajectory does not match topology!");
 146   out=ftp2FILE(efOUT,NFILE,fnm,"w");
 147   xvgr_legend(out,asize(leg),leg);
 148   nt=t;
 149   do {
 150     if (t >= nt) {
 151       rm_pbc(&(top->idef),top->atoms.nr,box,x,x_s);
 152       fp=ffopen(pdbfile,"w");
 153       write_pdbfile_indexed(fp,"Generated by do_shift",
 154                             atoms,x_s,box,' ',FALSE,gnx,index);
 155       fclose(fp);
 156
 157       if ((tot=popen(total,"w")) == NULL)
 158         perror("opening pipe to total");
 159       fprintf(tot,"%s\n",pdbfile);
 160       fprintf(tot,"%s\n",tmpfile);
 161       fprintf(tot,"3\n");
 162       fprintf(tot,"N\n");
 163       fprintf(tot,"%s\n",randf);
 164       fprintf(tot,"N\n");
 165       fprintf(tot,"N\n");
 166       if (pclose(tot) != 0)
 167         perror("closing pipe to total");
 168       cat(out,tmpfile,t);
 169       remove(pdbfile);
 170       remove(tmpfile);
 171       nt+=dt;
 172       nframe++;
 173     }
 174   } while(read_next_x(status,&t,natoms,x,box));
 175   close_trj(status);
 176   fclose(out);
 177
 178   thanx(stdout);
 179
 180   return 0;
 181 }