2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2006, The GROMACS development team,
6 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
54 #include "gmx_fatal.h"
60 void cat(FILE *out,char *fn,real t)
66 double f1,f2,f3,f4,f5,f6;
69 while ((ptr=fgets2(buf,255,in)) != NULL) {
70 sscanf(buf,"%d%d%s%s%lf%lf%lf%lf%lf%lf",
71 &anr,&rnr,rnm,anm,&f1,&f2,&f3,&f4,&f5,&f6);
72 fprintf(out,"%8g %10g %10g %10g %10g %10g %10g %s%d-%s%d\n",
73 t,f6,f1,f2,f3,f4,f5,rnm,rnr,anm,anr);
75 /*if ((int)strlen(buf) > 0)
76 fprintf(out,"%s\n",buf);*/
81 int main(int argc,char *argv[])
83 static char *desc[] = {
84 "[TT]do_shift[tt] reads a trajectory file and computes the chemical",
85 "shift for each time frame (or every [BB]dt[bb] ps) by",
86 "calling the 'total' program. If you do not have the total program,",
87 "get it. do_shift assumes that the total executable is in",
88 "[TT]/home/mdgroup/total/total[tt]. If that is not the case, then you should",
89 "set an environment variable [BB]TOTAL[bb] as in: [PAR]",
90 "[TT]setenv TOTAL /usr/local/bin/total[tt][PAR]",
91 "where the right hand side should point to the total executable.[PAR]",
92 "Output is printed in files [TT]shift.out[tt] where t is the time of the frame.[PAR]",
93 "The program also needs an input file called [BB]random.dat[bb] which",
94 "contains the random coil chemical shifts of all protons."
98 { "-dt", FALSE, etREAL, { &dt }, "Time interval between frames." }
100 static char *bugs[] = {
101 "The program is very slow"
103 static char *OXYGEN="O";
109 int i,natoms,nframe=0;
115 char pdbfile[32],tmpfile[32];
116 char total[256],*dptr;
118 { efTRX, "-f", NULL, ffREAD },
119 { efTPX, NULL, NULL, ffREAD },
120 { efNDX, NULL, NULL, ffREAD },
121 { efOUT, "-o", "shift", ffWRITE },
122 { efDAT, "-d", "random", ffREAD }
124 char *leg[] = { "shift","ring","anisCO","anisCN","sigmaE","sum" };
125 #define NFILE asize(fnm)
127 CopyRight(stdout,argv[0]);
128 parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_BE_NICE ,NFILE,fnm,
129 asize(pa),pa,asize(desc),desc,asize(bugs),bugs);
131 top=read_top(ftp2fn(efTPX,NFILE,fnm));
134 for(i=0; (i<atoms->nr); i++)
135 if ((strcmp(*atoms->atomname[i],"O1") == 0) ||
136 (strcmp(*atoms->atomname[i],"O2") == 0) ||
137 (strcmp(*atoms->atomname[i],"OXT") == 0) ||
138 (strcmp(*atoms->atomname[i],"OT") == 0))
139 atoms->atomname[i]=&OXYGEN;
140 rd_index(ftp2fn(efNDX,NFILE,fnm),1,&gnx,&index,&grpnm);
144 strcpy(pdbfile,"dsXXXXXX");
146 strcpy(tmpfile,"dsXXXXXX");
148 fprintf(stderr,"pdbfile = %s\ntmpfile = %s\n",pdbfile,tmpfile);
150 if ((dptr=getenv("TOTAL")) == NULL)
151 dptr="/home/mdgroup/total/total";
152 sprintf(total,"%s > /dev/null",dptr);
153 fprintf(stderr,"total cmd='%s'\n",total);
154 randf=ftp2fn(efDAT,NFILE,fnm);
156 natoms=read_first_x(&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
157 if (natoms != atoms->nr)
158 gmx_fatal(FARGS,"Trajectory does not match topology!");
159 out=ftp2FILE(efOUT,NFILE,fnm,"w");
160 xvgr_legend(out,asize(leg),leg);
164 rm_pbc(&(top->idef),top->atoms.nr,box,x,x_s);
165 fp=ffopen(pdbfile,"w");
166 write_pdbfile_indexed(fp,"Generated by do_shift",
167 atoms,x_s,box,0,-1,gnx,index);
170 if ((tot=popen(total,"w")) == NULL)
171 perror("opening pipe to total");
172 fprintf(tot,"%s\n",pdbfile);
173 fprintf(tot,"%s\n",tmpfile);
176 fprintf(tot,"%s\n",randf);
179 if (pclose(tot) != 0)
180 perror("closing pipe to total");
187 } while(read_next_x(status,&t,natoms,x,box));