src/contrib/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 3.0
  11  *
  12  * Copyright (c) 1991-2001
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * Do check out http://www.gromacs.org , or mail us at gromacs@gromacs.org .
  32  *
  33  * And Hey:
  34  * Gromacs Runs On Most of All Computer Systems
  35  */
  36 static char *SRCID_do_shift_c = "$Id$";
  37 #include <stdlib.h>
  38 #include "errno.h"
  39 #include "sysstuff.h"
  40 #include "typedefs.h"
  41 #include "string2.h"
  42 #include "strdb.h"
  43 #include "macros.h"
  44 #include "smalloc.h"
  45 #include "mshift.h"
  46 #include "statutil.h"
  47 #include "copyrite.h"
  48 #include "confio.h"
  49 #include "fatal.h"
  50 #include "xvgr.h"
  51 #include "gstat.h"
  52 #include "rdgroup.h"
  53 #include "pdbio.h"
  54
  55 void cat(FILE *out,char *fn,real t)
  56 {
  57   FILE *in;
  58   char *ptr,buf[256];
  59   int    anr,rnr;
  60   char   anm[24],rnm[24];
  61   double f1,f2,f3,f4,f5,f6;
  62
  63   in=ffopen(fn,"r");
  64   while ((ptr=fgets2(buf,255,in)) != NULL) {
  65     sscanf(buf,"%d%d%s%s%lf%lf%lf%lf%lf%lf",
  66            &anr,&rnr,rnm,anm,&f1,&f2,&f3,&f4,&f5,&f6);
  67     fprintf(out,"%8g  %10g  %10g  %10g  %10g  %10g  %10g  %s%d-%s%d\n",
  68             t,f6,f1,f2,f3,f4,f5,rnm,rnr,anm,anr);
  69   }
  70   /*if ((int)strlen(buf) > 0)
  71     fprintf(out,"%s\n",buf);*/
  72   fflush(out);
  73   fclose(in);
  74 }
  75
  76 int main(int argc,char *argv[])
  77 {
  78   static char *desc[] = {
  79     "do_shift reads a trajectory file and computes the chemical",
  80     "shift for each time frame (or every [BB]dt[bb] ps) by",
  81     "calling the 'total' program. If you do not have the total program,",
  82     "get it. do_shift assumes that the total executable is in",
  83     "/home/mdgroup/total/total. If that is not the case, then you should",
  84     "set an environment variable [BB]TOTAL[bb] as in: [PAR]",
  85     "[TT]setenv TOTAL /usr/local/bin/total[tt][PAR]",
  86     "where the right hand side should point to the total executable.[PAR]",
  87     "Output is printed in files shift.out where t is the time of the frame.[PAR]",
  88     "The program also needs an input file called [BB]random.dat[bb] which",
  89     "contains the random coil chemical shifts of all protons."
  90   };
  91   static real dt=0.0;
  92   t_pargs pa[] = {
  93     { "-dt", FALSE, etREAL, &dt, "Time interval between frames." }
  94   };
  95   static char *bugs[] = {
  96     "The program is very slow"
  97   };
  98   static     char *OXYGEN="O";
  99   FILE       *out,*tot,*fp;
 100   t_topology *top;
 101   t_atoms    *atoms;
 102   int        status,nres;
 103   real       t,nt;
 104   int        i,natoms,nframe=0;
 105   matrix     box;
 106   int        gnx;
 107   char       *grpnm,*randf;
 108   atom_id    *index;
 109   rvec       *x,*x_s;
 110   char       pdbfile[L_tmpnam],tmpfile[L_tmpnam];
 111   char       total[256],*dptr;
 112   t_filenm   fnm[] = {
 113     { efTRX, "-f",   NULL,     ffREAD },
 114     { efTPX, NULL,   NULL,     ffREAD },
 115     { efNDX, NULL,   NULL,     ffREAD },
 116     { efOUT, "-o",   "shift",  ffWRITE },
 117     { efDAT, "-d",   "random", ffREAD }
 118   };
 119   char *leg[] = { "shift","ring","anisCO","anisCN","sigmaE","sum" };
 120 #define NFILE asize(fnm)
 121
 122   CopyRight(stdout,argv[0]);
 123   parse_common_args(&argc,argv,PCA_CAN_TIME,TRUE,NFILE,fnm,
 124                     asize(pa),pa,asize(desc),desc,asize(bugs),bugs);
 125
 126   top=read_top(ftp2fn(efTPX,NFILE,fnm));
 127   atoms=&(top->atoms);
 128   nres=atoms->nres;
 129   for(i=0; (i<atoms->nr); i++)
 130     if ((strcmp(*atoms->atomname[i],"O1") == 0) ||
 131         (strcmp(*atoms->atomname[i],"O2") == 0) ||
 132         (strcmp(*atoms->atomname[i],"OXT") == 0) ||
 133         (strcmp(*atoms->atomname[i],"OT") == 0))
 134       atoms->atomname[i]=&OXYGEN;
 135   rd_index(ftp2fn(efNDX,NFILE,fnm),1,&gnx,&index,&grpnm);
 136
 137   snew(x_s,atoms->nr);
 138
 139   (void) tmpnam(pdbfile);
 140   (void) tmpnam(tmpfile);
 141   fprintf(stderr,"pdbfile = %s\ntmpfile = %s\n",pdbfile,tmpfile);
 142
 143   if ((dptr=getenv("TOTAL")) == NULL)
 144     dptr="/home/mdgroup/total/total";
 145   sprintf(total,"%s > /dev/null",dptr);
 146   fprintf(stderr,"total cmd='%s'\n",total);
 147   randf=ftp2fn(efDAT,NFILE,fnm);
 148
 149   natoms=read_first_x(&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
 150   if (natoms != atoms->nr)
 151     fatal_error(0,"Trajectory does not match topology!");
 152   out=ftp2FILE(efOUT,NFILE,fnm,"w");
 153   xvgr_legend(out,asize(leg),leg);
 154   nt=t;
 155   do {
 156     if (t >= nt) {
 157       rm_pbc(&(top->idef),top->atoms.nr,box,x,x_s);
 158       fp=ffopen(pdbfile,"w");
 159       write_pdbfile_indexed(fp,"Generated by do_shift",
 160                             atoms,x_s,box,0,-1,gnx,index);
 161       fclose(fp);
 162
 163       if ((tot=popen(total,"w")) == NULL)
 164         perror("opening pipe to total");
 165       fprintf(tot,"%s\n",pdbfile);
 166       fprintf(tot,"%s\n",tmpfile);
 167       fprintf(tot,"3\n");
 168       fprintf(tot,"N\n");
 169       fprintf(tot,"%s\n",randf);
 170       fprintf(tot,"N\n");
 171       fprintf(tot,"N\n");
 172       if (pclose(tot) != 0)
 173         perror("closing pipe to total");
 174       cat(out,tmpfile,t);
 175       remove(pdbfile);
 176       remove(tmpfile);
 177       nt+=dt;
 178       nframe++;
 179     }
 180   } while(read_next_x(status,&t,natoms,x,box));
 181   close_trj(status);
 182   fclose(out);
 183
 184   thanx(stderr);
 185
 186   return 0;
 187 }