1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
3 * This source code is NOT REALLY part of
7 * GROningen MAchine for Chemical Simulations
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-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
32 * Author: do_multiprot was written by Ran Friedman <r.friedman@bioc.uzh.ch>
35 * Green Red Orange Magenta Azure Cyan Skyblue
44 #include "gromacs/utility/smalloc.h"
46 #include "gromacs/commandline/pargs.h"
48 #include "gromacs/fileio/pdbio.h"
49 #include "gromacs/utility/fatalerror.h"
50 #include "gromacs/fileio/xvgr.h"
51 #include "gromacs/fileio/matio.h"
54 #include "gromacs/fileio/tpxio.h"
58 #include "gromacs/fileio/confio.h"
59 #include "gromacs/fileio/gmxfio.h"
66 static void process_multiprot_output(const char *fn, real *rmsd, int *nres, rvec rotangles,
67 rvec translation, bool bCountres, t_countres *countres)
76 mpoutput=gmx_ffopen (fn,"r");
80 fgets(line, 256, mpoutput);
81 } while (strstr(line,"Match List") == NULL);
82 fgets(line, 256, mpoutput);
84 string = strtok (line,".");
85 while (i<2 && string != NULL) {
86 string = strtok (NULL,". ");
92 while (countres[j].resnr!=res)
96 fgets(line, 256, mpoutput);
97 } while (strstr(line,"End of Match List") == NULL);
101 fgets(line, 256, mpoutput);
102 } while (strstr(line,"Trans : ") == NULL);
104 string = strtok (line," :");
105 string = strtok (NULL," ");
106 while (i<3 && string != NULL) {
107 string = strtok (NULL," ");
108 rotangles[i]=atof(string);
112 while (i<3 && string != NULL) {
113 string = strtok (NULL," ");
114 translation[i]=atof(string)/10;
118 gmx_warning("Not enough values for rotation and translation vectors in the output of multiprot");
121 rotangles[YY]=rotangles[YY]*(-1);
124 fgets(line, 256, mpoutput);
125 if (strstr(line,"RMSD : ") != NULL) {
126 string = strtok (line,":");
127 string = strtok (NULL,":");
128 (*rmsd)=atof(string)/10;
132 fgets(line,256, mpoutput);
133 if (strstr(line,"Match List") != NULL) {
134 string = strtok (line,":");
135 string = strtok (NULL,":");
136 (*nres) = atoi(string);
139 gmx_ffclose(mpoutput);
142 int main(int argc,char *argv[])
144 const char *desc[] = {
145 "[TT]do_multiprot[tt] ",
146 "reads a trajectory file and aligns it to a reference structure ",
148 "calling the multiprot program. This allows you to use a reference",
149 "structure whose sequence is different than that of the protein in the ",
150 "trajectory, since the alignment is based on the geometry, not sequence.",
151 "The output of [TT]do_multiprot[tt] includes the rmsd and the number of residues",
152 "on which it was calculated.",
154 "An aligned trajectory file is generated with the [TT]-ox[tt] option.[PAR]",
155 "With the [TT]-cr[tt] option, the number of hits in the alignment is given",
156 "per residue. This number can be between 0 and the number of frames, and",
157 "indicates the structural conservation of this residue.[PAR]",
158 "If you do not have the [TT]multiprot[tt] program, get it. [TT]do_multiprot[tt] assumes",
159 "that the [TT]multiprot[tt] executable is [TT]/usr/local/bin/multiprot[tt]. If this is ",
160 "not the case, then you should set an environment variable [BB]MULTIPROT[bb]",
161 "pointing to the [TT]multiprot[tt] executable, e.g.: [PAR]",
162 "[TT]setenv MULTIPROT /usr/MultiProtInstall/multiprot.Linux[tt][PAR]",
163 "Note that at the current implementation only binary alignment (your",
164 "molecule to a reference) is supported. In addition, note that the ",
165 "by default [TT]multiprot[tt] aligns the two proteins on their C-alpha carbons.",
166 "and that this depends on the [TT]multiprot[tt] parameters which are not dealt ",
167 "with here. Thus, the C-alpha carbons is expected to give the same "
168 "results as choosing the whole protein and will be slightly faster.[PAR]",
169 "For information about [TT]multiprot[tt], see:",
170 "http://bioinfo3d.cs.tau.ac.il/MultiProt/.[PAR]"
172 static bool bVerbose;
174 { "-v", FALSE, etBOOL, {&bVerbose},
175 "HIDDENGenerate miles of useless information" }
178 const char *bugs[] = {
179 "The program is very slow, since multiprot is run externally"
183 t_trxstatus *trxout=NULL;
184 FILE *tapein,*fo,*frc,*tmpf,*out=NULL,*fres=NULL;
186 const char *fn="2_sol.res";
189 t_atoms *atoms,ratoms,useatoms;
194 int i,j,natoms,nratoms,nframe=0,model_nr=-1;
195 int cur_res,prev_res;
197 t_countres *countres=NULL;
203 char pdbfile[32],refpdb[256],title[256],rtitle[256],filemode[5];
205 char multiprot[256],*mptr;
209 bool bTrjout,bCountres;
210 const char *TrjoutFile=NULL;
212 static rvec translation={0,0,0},rotangles={0,0,0};
213 gmx_rmpbc_t gpbc=NULL;
216 { efTRX, "-f", NULL, ffREAD },
217 { efTPS, NULL, NULL, ffREAD },
218 { efNDX, NULL, NULL, ffOPTRD },
219 { efSTX, "-r", NULL , ffREAD },
220 { efXVG, "-o", "rmss", ffWRITE },
221 { efXVG, "-rc", "rescount", ffWRITE},
222 { efXVG, "-cr", "countres", ffOPTWR},
223 { efTRX, "-ox", "aligned", ffOPTWR }
225 #define NFILE asize(fnm)
227 CopyRight(stderr,argv[0]);
228 parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_CAN_VIEW | PCA_TIME_UNIT | PCA_BE_NICE ,
229 NFILE,fnm, asize(pa),pa, asize(desc),desc,
230 asize(bugs),bugs,&oenv
232 fnRef=opt2fn("-r",NFILE,fnm);
233 bTrjout = opt2bSet("-ox",NFILE,fnm);
234 bCountres= opt2bSet("-cr",NFILE,fnm);
237 TrjoutFile = opt2fn_null("-ox",NFILE,fnm);
240 read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,&xp,NULL,box,FALSE);
241 gpbc = gmx_rmpbc_init(&top.idef,ePBC,top.atoms.nr,box);
246 get_stx_coordnum(fnRef,&nratoms);
247 init_t_atoms(&ratoms,nratoms,TRUE);
249 read_stx_conf(fnRef,rtitle,&ratoms,xr,NULL,&ePBC,rbox);
252 fprintf(stderr,"Read %d atoms\n",atoms->nr);
253 fprintf(stderr,"Read %d reference atoms\n",ratoms.nr);
256 snew(countres,ratoms.nres);
259 for (i=0;i<ratoms.nr;i++) {
261 cur_res=ratoms.atom[i].resind;
262 if (cur_res != prev_res) {
263 countres[j].resnr=cur_res;
269 get_index(atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&gnx,&index,&grpnm);
272 for(i=0; (i<gnx); i++) {
273 if (atoms->atom[index[i]].resind != nr0) {
274 nr0=atoms->atom[index[i]].resind;
278 fprintf(stderr,"There are %d residues in your selected group\n",nres);
280 strcpy(pdbfile,"ddXXXXXX");
282 if ((tmpf = fopen(pdbfile,"w")) == NULL) {
283 sprintf(pdbfile,"%ctmp%cfilterXXXXXX",DIR_SEPARATOR,DIR_SEPARATOR);
285 if ((tmpf = fopen(pdbfile,"w")) == NULL) {
286 gmx_fatal(FARGS,"Can not open tmp file %s",pdbfile);
294 strcpy(refpdb,"ddXXXXXX");
296 strcat(refpdb,".pdb");
297 write_sto_conf(refpdb,rtitle,&ratoms,xr,NULL,ePBC,rbox);
300 strcpy(refpdb,fnRef);
303 if ((mptr=getenv("MULTIPROT")) == NULL) {
304 mptr="/usr/local/bin/multiprot";
306 if (!gmx_fexist(mptr)) {
307 gmx_fatal(FARGS,"MULTIPROT executable (%s) does not exist (use setenv MULTIPROT)",
310 sprintf (multiprot,"%s %s %s > /dev/null %s",
311 mptr, refpdb, pdbfile, "2> /dev/null");
314 fprintf(stderr,"multiprot cmd='%s'\n",multiprot);
316 if (!read_first_frame(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&fr,TRX_READ_X))
317 gmx_fatal(FARGS,"Could not read a frame from %s",ftp2fn(efTRX,NFILE,fnm));
323 outftp=fn2ftp(TrjoutFile);
325 fprintf(stderr,"Will write %s: %s\n",ftp2ext(ftp),ftp2desc(outftp));
326 strcpy(filemode,"w");
333 trxout = open_trx(TrjoutFile,filemode);
338 /* Make atoms struct for output in GRO or PDB files */
339 /* get memory for stuff to go in pdb file */
340 init_t_atoms(&useatoms,nout,FALSE);
341 sfree(useatoms.resinfo);
342 useatoms.resinfo=atoms->resinfo;
343 for(i=0;(i<nout);i++) {
344 useatoms.atomname[i]=atoms->atomname[i];
345 useatoms.atom[i]=atoms->atom[i];
346 useatoms.nres=max(useatoms.nres,useatoms.atom[i].resind+1);
349 out=gmx_ffopen(TrjoutFile,filemode);
353 fprintf(gmx_fio_getfp(trx_get_fileio(trxout)),"Generated by %s. #atoms=%d, a BOX is"
354 " stored in this file.\n",ShortProgram(),nout);
357 if (natoms > atoms->nr) {
358 gmx_fatal(FARGS,"\nTrajectory does not match topology!");
361 gmx_fatal(FARGS,"\nTrajectory does not match selected group!");
364 fo = xvgropen(opt2fn("-o",NFILE,fnm),"RMSD","Time (ps)","RMSD (nm)",oenv);
365 frc = xvgropen(opt2fn("-rc",NFILE,fnm),"Number of Residues in the alignment","Time (ps)","Residues",oenv);
368 t = output_env_conv_time(oenv,fr.time);
369 gmx_rmpbc(gpbc,natoms,fr.box,fr.x);
370 tapein=gmx_ffopen(pdbfile,"w");
371 write_pdbfile_indexed(tapein,NULL,atoms,fr.x,ePBC,fr.box,' ',-1,gnx,index,NULL,TRUE);
375 process_multiprot_output(fn, &rmsd, &nres2,rotangles,translation,bCountres,countres);
376 fprintf(fo,"%12.7f",t);
377 fprintf(fo," %12.7f\n",rmsd);
378 fprintf(frc,"%12.7f",t);
379 fprintf(frc,"%12d\n",nres2);
381 rotate_conf(natoms,fr.x,NULL,rotangles[XX],rotangles[YY],rotangles[ZZ]);
382 for(i=0; i<natoms; i++) {
383 rvec_inc(fr.x[i],translation);
390 write_trxframe(trxout,&fr,NULL);
395 sprintf(out_title,"Generated by do_multiprot : %s t= %g %s",
396 title,output_env_conv_time(oenv,fr.time),output_env_get_time_unit(oenv));
399 write_hconf_p(out,out_title,&useatoms,prec2ndec(fr.prec),
403 fprintf(out,"REMARK GENERATED BY DO_MULTIPROT\n");
404 sprintf(out_title,"%s t= %g %s",title,output_env_conv_time(oenv,fr.time),output_env_get_time_unit(oenv));
405 /* if reading from pdb, we want to keep the original
406 model numbering else we write the output frame
407 number plus one, because model 0 is not allowed in pdb */
408 if (ftp==efPDB && fr.step > model_nr) {
414 write_pdbfile(out,out_title,&useatoms,fr.x,ePBC,fr.box,' ',model_nr,NULL,TRUE);
417 fr.title = out_title;
418 fr.bTitle = (nframe == 0);
422 write_g96_conf(out,&fr,-1,NULL);
428 } while(read_next_frame(oenv,status,&fr));
430 fres= xvgropen(opt2fn("-cr",NFILE,fnm),"Number of frames in which the residues are aligned to","Residue","Number",oenv);
431 for (i=0;i<ratoms.nres;i++) {
432 fprintf(fres,"%10d %12d\n",countres[i].resnr,countres[i].count);
438 fprintf(stderr,"\n");
440 if (trxout != NULL) {
443 else if (out != NULL) {
446 view_all(oenv,NFILE, fnm);
451 free_t_atoms(&ratoms,TRUE);
453 if (outftp==efPDB || outftp==efGRO || outftp==efG96) {
454 free_t_atoms(&useatoms,TRUE);