383a0dea12ab51ef464aeaa6910acc36b94528e8
[alexxy/gromacs.git] / 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 "typedefs.h"
42 #include "gromacs/utility/cstringutil.h"
43 #include "macros.h"
44 #include "gromacs/utility/smalloc.h"
45 #include "mshift.h"
46 #include "gromacs/commandline/pargs.h"
47 #include "copyrite.h"
48 #include "gromacs/fileio/confio.h"
49 #include "gromacs/utility/fatalerror.h"
50 #include "xvgr.h"
51 #include "gstat.h"
52 #include "index.h"
53 #include "gromacs/fileio/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=gmx_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   gmx_ffclose(in);
74 }
75
76 int main(int argc,char *argv[])
77 {
78   static char *desc[] = {
79     "[TT]do_shift[tt] 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     "[TT]/home/mdgroup/total/total[tt]. 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 [TT]shift.out[tt] 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[32],tmpfile[32];
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 | PCA_BE_NICE ,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   strcpy(pdbfile,"dsXXXXXX");
140   gmx_tmpnam(pdbfile);
141   strcpy(tmpfile,"dsXXXXXX");
142   gmx_tmpnam(tmpfile);
143   fprintf(stderr,"pdbfile = %s\ntmpfile = %s\n",pdbfile,tmpfile);
144   
145   if ((dptr=getenv("TOTAL")) == NULL)
146     dptr="/home/mdgroup/total/total";
147   sprintf(total,"%s > /dev/null",dptr);
148   fprintf(stderr,"total cmd='%s'\n",total);
149   randf=ftp2fn(efDAT,NFILE,fnm);
150   
151   natoms=read_first_x(&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
152   if (natoms != atoms->nr) 
153     gmx_fatal(FARGS,"Trajectory does not match topology!");
154   out=ftp2FILE(efOUT,NFILE,fnm,"w");
155   xvgr_legend(out,asize(leg),leg);
156   nt=t;
157   do {
158     if (t >= nt) {
159       rm_pbc(&(top->idef),top->atoms.nr,box,x,x_s);
160       fp=gmx_ffopen(pdbfile,"w");
161       write_pdbfile_indexed(fp,"Generated by do_shift",
162                             atoms,x_s,box,0,-1,gnx,index);
163       gmx_ffclose(fp);
164       
165       if ((tot=popen(total,"w")) == NULL)
166         perror("opening pipe to total");
167       fprintf(tot,"%s\n",pdbfile);
168       fprintf(tot,"%s\n",tmpfile);
169       fprintf(tot,"3\n");
170       fprintf(tot,"N\n");
171       fprintf(tot,"%s\n",randf);
172       fprintf(tot,"N\n");
173       fprintf(tot,"N\n");
174       if (pclose(tot) != 0)
175         perror("closing pipe to total");
176       cat(out,tmpfile,t);
177       remove(pdbfile);
178       remove(tmpfile);
179       nt+=dt;
180       nframe++;
181     }
182   } while(read_next_x(status,&t,natoms,x,box));
183   close_trj(status);
184   gmx_ffclose(out);
185   
186   gmx_thanx(stderr);
187   
188   return 0;
189 }