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