4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
12 * Copyright (c) 1991-1999
13 * BIOSON Research Institute, Dept. of Biophysical Chemistry
14 * University of Groningen, The Netherlands
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)
21 * Also check out our WWW page:
22 * http://md.chem.rug.nl/~gmx
27 * Great Red Oystrich Makes All Chemists Sane
29 static char *SRCID_do_shift_c = "$Id$";
49 void cat(FILE *out,char *fn,real t)
55 double f1,f2,f3,f4,f5,f6;
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);
64 /*if ((int)strlen(buf) > 0)
65 fprintf(out,"%s\n",buf);*/
70 int main(int argc,char *argv[])
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."
87 { "-dt", FALSE, etREAL, &dt, "Time interval between frames." }
89 static char *bugs[] = {
90 "The program is very slow"
92 static char *OXYGEN="O";
98 int i,natoms,nframe=0;
104 char pdbfile[L_tmpnam],tmpfile[L_tmpnam];
105 char total[256],*dptr;
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 }
113 char *leg[] = { "shift","ring","anisCO","anisCN","sigmaE","sum" };
114 #define NFILE asize(fnm)
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);
120 top=read_top(ftp2fn(efTPX,NFILE,fnm));
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);
133 (void) tmpnam(pdbfile);
134 (void) tmpnam(tmpfile);
135 fprintf(stderr,"pdbfile = %s\ntmpfile = %s\n",pdbfile,tmpfile);
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);
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);
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,0,-1,gnx,index);
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);
163 fprintf(tot,"%s\n",randf);
166 if (pclose(tot) != 0)
167 perror("closing pipe to total");
174 } while(read_next_x(status,&t,natoms,x,box));