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-2004, 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.
64 static void calc_dist(int nind,atom_id index[],rvec x[],int ePBC,matrix box,
73 set_pbc(&pbc,ePBC,box);
74 for(i=0; (i<nind-1); i++) {
76 for(j=i+1; (j<nind); j++) {
77 pbc_dx(&pbc,xi,x[index[j]],dx);
84 static void calc_dist_tot(int nind,atom_id index[],rvec x[],
86 real **d, real **dtot, real **dtot2,
87 gmx_bool bNMR, real **dtot1_3, real **dtot1_6)
91 real temp, temp2, temp1_3;
95 set_pbc(&pbc,ePBC,box);
96 for(i=0; (i<nind-1); i++) {
98 for(j=i+1; (j<nind); j++) {
99 pbc_dx(&pbc,xi,x[index[j]],dx);
106 temp1_3 = 1.0/(temp*temp2);
107 dtot1_3[i][j] += temp1_3;
108 dtot1_6[i][j] += temp1_3*temp1_3;
114 static void calc_nmr(int nind, int nframes, real **dtot1_3, real **dtot1_6,
115 real *max1_3, real *max1_6)
118 real temp1_3,temp1_6;
120 for(i=0; (i<nind-1); i++) {
121 for(j=i+1; (j<nind); j++) {
122 temp1_3 = pow(dtot1_3[i][j]/nframes,-1.0/3.0);
123 temp1_6 = pow(dtot1_6[i][j]/nframes,-1.0/6.0);
124 if (temp1_3 > *max1_3) *max1_3 = temp1_3;
125 if (temp1_6 > *max1_6) *max1_6 = temp1_6;
126 dtot1_3[i][j] = temp1_3;
127 dtot1_6[i][j] = temp1_6;
128 dtot1_3[j][i] = temp1_3;
129 dtot1_6[j][i] = temp1_6;
134 static char Hnum[] = "123";
159 static int read_equiv(const char *eq_fn, t_equiv ***equivptr)
162 char line[STRLEN],resname[10],atomname[10],*lp;
163 int neq, na, n, resnr;
166 fp = ffopen(eq_fn,"r");
169 while(get_a_line(fp,line,STRLEN)) {
171 /* this is not efficient, but I'm lazy */
175 if (sscanf(lp,"%s %n",atomname,&n)==1) {
178 equiv[neq][0].nname=strdup(atomname);
179 while (sscanf(lp,"%d %s %s %n",&resnr,resname,atomname,&n)==3) {
180 /* this is not efficient, but I'm lazy (again) */
181 srenew(equiv[neq], na+1);
182 equiv[neq][na].rnr=resnr-1;
183 equiv[neq][na].rname=strdup(resname);
184 equiv[neq][na].aname=strdup(atomname);
185 if (na>0) equiv[neq][na].nname=NULL;
190 /* make empty element as flag for end of array */
191 srenew(equiv[neq], na+1);
192 equiv[neq][na].rnr=NOTSET;
193 equiv[neq][na].rname=NULL;
194 equiv[neq][na].aname=NULL;
206 static void dump_equiv(FILE *out, int neq, t_equiv **equiv)
210 fprintf(out,"Dumping equivalent list\n");
211 for (i=0; i<neq; i++) {
212 fprintf(out,"%s",equiv[i][0].nname);
213 for(j=0; equiv[i][j].rnr!=NOTSET; j++)
214 fprintf(out," %d %s %s",
215 equiv[i][j].rnr,equiv[i][j].rname,equiv[i][j].aname);
220 static gmx_bool is_equiv(int neq, t_equiv **equiv, char **nname,
221 int rnr1, char *rname1, char *aname1,
222 int rnr2, char *rname2, char *aname2)
228 /* we can terminate each loop when bFound is true! */
229 for (i=0; i<neq && !bFound; i++) {
230 /* find first atom */
231 for(j=0; equiv[i][j].rnr!=NOTSET && !bFound; j++)
232 bFound = ( equiv[i][j].rnr==rnr1 &&
233 strcmp(equiv[i][j].rname, rname1)==0 &&
234 strcmp(equiv[i][j].aname, aname1)==0 );
236 /* find second atom */
238 for(j=0; equiv[i][j].rnr!=NOTSET && !bFound; j++)
239 bFound = ( equiv[i][j].rnr==rnr2 &&
240 strcmp(equiv[i][j].rname, rname2)==0 &&
241 strcmp(equiv[i][j].aname, aname2)==0 );
245 *nname = strdup(equiv[i-1][0].nname);
250 static int analyze_noe_equivalent(const char *eq_fn,
251 t_atoms *atoms, int isize, atom_id *index,
253 atom_id *noe_index, t_noe_gr *noe_gr)
256 int i, j, anmil, anmjl, rnri, rnrj, gi, groupnr, neq;
257 char *anmi, *anmj, **nnm;
258 gmx_bool bMatch,bEquiv;
264 neq=read_equiv(eq_fn,&equiv);
265 if (debug) dump_equiv(debug,neq,equiv);
272 for(i=0; i<isize; i++) {
273 if (equiv && i<isize-1) {
274 /* check explicit list of equivalent atoms */
277 rnri=atoms->atom[index[i]].resind;
278 rnrj=atoms->atom[index[j]].resind;
280 is_equiv(neq, equiv, &nnm[i],
281 rnri,*atoms->resinfo[rnri].name,*atoms->atomname[index[i]],
282 rnrj,*atoms->resinfo[rnrj].name,*atoms->atomname[index[j]]);
284 nnm[j]=strdup(nnm[i]);
286 /* set index for matching atom */
287 noe_index[j]=groupnr;
288 /* skip matching atom */
291 } while ( bEquiv && i<isize-1 );
295 /* look for triplets of consecutive atoms with name XX?,
296 X are any number of letters or digits and ? goes from 1 to 3
297 This is supposed to cover all CH3 groups and the like */
298 anmi = *atoms->atomname[index[i]];
299 anmil = strlen(anmi);
300 bMatch = i<isize-3 && anmi[anmil-1]=='1';
303 anmj = *atoms->atomname[index[i+j]];
304 anmjl = strlen(anmj);
305 bMatch = bMatch && ( anmil==anmjl && anmj[anmjl-1]==Hnum[j] &&
306 strncmp(anmi, anmj, anmil-1)==0 );
308 /* set index for this atom */
309 noe_index[i]=groupnr;
311 /* set index for next two matching atoms */
313 noe_index[i+j]=groupnr;
314 /* skip matching atoms */
321 /* make index without looking for equivalent atoms */
322 for(i=0; i<isize; i++)
326 noe_index[isize]=groupnr;
330 for(i=0; i<isize; i++) {
331 rnri=atoms->atom[index[i]].resind;
332 fprintf(debug,"%s %s %d -> %s\n",*atoms->atomname[index[i]],
333 *atoms->resinfo[rnri].name,rnri,nnm[i]?nnm[i]:"");
336 for(i=0; i<isize; i++) {
338 if (!noe_gr[gi].aname) {
340 noe_gr[gi].anr=index[i];
342 noe_gr[gi].aname=strdup(nnm[i]);
344 noe_gr[gi].aname=strdup(*atoms->atomname[index[i]]);
345 if ( noe_index[i]==noe_index[i+1] )
346 noe_gr[gi].aname[strlen(noe_gr[gi].aname)-1]='*';
348 noe_gr[gi].rnr=atoms->atom[index[i]].resind;
349 noe_gr[gi].rname=strdup(*atoms->resinfo[noe_gr[gi].rnr].name);
350 /* dump group definitions */
351 if (debug) fprintf(debug,"%d %d %d %d %s %s %d\n",i,gi,
352 noe_gr[gi].ianr,noe_gr[gi].anr,noe_gr[gi].aname,
353 noe_gr[gi].rname,noe_gr[gi].rnr);
356 for(i=0; i<isize; i++)
363 /* #define NSCALE 3 */
364 /* static char *noe_scale[NSCALE+1] = { */
365 /* "strong", "medium", "weak", "none" */
369 static char *noe2scale(real r3, real r6, real rmax)
372 static char buf[NSCALE+1];
374 /* r goes from 0 to rmax
375 NSCALE*r/rmax goes from 0 to NSCALE
376 NSCALE - NSCALE*r/rmax goes from NSCALE to 0 */
377 s3 = NSCALE - min(NSCALE, (int)(NSCALE*r3/rmax));
378 s6 = NSCALE - min(NSCALE, (int)(NSCALE*r6/rmax));
389 static void calc_noe(int isize, atom_id *noe_index,
390 real **dtot1_3, real **dtot1_6, int gnr, t_noe **noe)
394 /* make half matrix of noe-group distances from atom distances */
395 for(i=0; i<isize; i++) {
397 for(j=i; j<isize; j++) {
400 noe[gi][gj].i_3+=pow(dtot1_3[i][j],-3);
401 noe[gi][gj].i_6+=pow(dtot1_6[i][j],-6);
407 for(j=i+1; j<gnr; j++) {
408 noe[i][j].r_3 = pow(noe[i][j].i_3/noe[i][j].nr,-1.0/3.0);
409 noe[i][j].r_6 = pow(noe[i][j].i_6/noe[i][j].nr,-1.0/6.0);
410 noe[j][i] = noe[i][j];
414 static void write_noe(FILE *fp,int gnr,t_noe **noe,t_noe_gr *noe_gr,real rmax)
417 real r3, r6, min3, min6;
418 char buf[10],b3[10],b6[10];
423 ";%4s %3s %4s %4s%3s %4s %4s %4s %4s%3s %5s %5s %8s %2s %2s %s\n",
424 "ianr","anr","anm","rnm","rnr","ianr","anr","anm","rnm","rnr",
425 "1/r^3","1/r^6","intnsty","Dr","Da","scale");
426 for(i=0; i<gnr; i++) {
428 for(j=i+1; j<gnr; j++) {
434 if ( r3 < rmax || r6 < rmax ) {
435 if (grj.rnr == gri.rnr)
436 sprintf(buf,"%2d", grj.anr-gri.anr);
440 sprintf(b3,"%-5.3f",r3);
444 sprintf(b6,"%-5.3f",r6);
448 "%4d %4d %4s %4s%3d %4d %4d %4s %4s%3d %5s %5s %8d %2d %2s %s\n",
449 gri.ianr+1, gri.anr+1, gri.aname, gri.rname, gri.rnr+1,
450 grj.ianr+1, grj.anr+1, grj.aname, grj.rname, grj.rnr+1,
451 b3, b6, (int)(noe[i][j].i_6+0.5), grj.rnr-gri.rnr, buf,
452 noe2scale(r3, r6, rmax));
456 #define MINI ((i==3)?min3:min6)
459 fprintf(stdout,"NOTE: no 1/r^%d averaged distances found below %g, "
460 "smallest was %g\n", i, rmax, MINI );
462 fprintf(stdout,"Smallest 1/r^%d averaged distance was %g\n", i, MINI );
466 static void calc_rms(int nind, int nframes,
467 real **dtot, real **dtot2,
468 real **rmsmat, real *rmsmax,
469 real **rmscmat, real *rmscmax,
470 real **meanmat, real *meanmax)
473 real mean, mean2, rms, rmsc;
474 /* N.B. dtot and dtot2 contain the total distance and the total squared
475 * distance respectively, BUT they return RMS and the scaled RMS resp.
482 for(i=0; (i<nind-1); i++) {
483 for(j=i+1; (j<nind); j++) {
484 mean =dtot[i][j]/nframes;
485 mean2=dtot2[i][j]/nframes;
486 rms=sqrt(max(0,mean2-mean*mean));
488 if (mean > *meanmax) *meanmax=mean;
489 if (rms > *rmsmax ) *rmsmax =rms;
490 if (rmsc > *rmscmax) *rmscmax=rmsc;
491 meanmat[i][j]=meanmat[j][i]=mean;
492 rmsmat[i][j] =rmsmat[j][i] =rms;
493 rmscmat[i][j]=rmscmat[j][i]=rmsc;
498 real rms_diff(int natom,real **d,real **d_r)
504 for(i=0; (i<natom-1); i++)
505 for(j=i+1; (j<natom); j++) {
509 r2/=(natom*(natom-1))/2;
514 int gmx_rmsdist (int argc,char *argv[])
516 const char *desc[] = {
517 "[TT]g_rmsdist[tt] computes the root mean square deviation of atom distances,",
518 "which has the advantage that no fit is needed like in standard RMS",
519 "deviation as computed by [TT]g_rms[tt].",
520 "The reference structure is taken from the structure file.",
521 "The RMSD at time t is calculated as the RMS",
522 "of the differences in distance between atom-pairs in the reference",
523 "structure and the structure at time t.[PAR]",
524 "[TT]g_rmsdist[tt] can also produce matrices of the rms distances, rms distances",
525 "scaled with the mean distance and the mean distances and matrices with",
526 "NMR averaged distances (1/r^3 and 1/r^6 averaging). Finally, lists",
527 "of atom pairs with 1/r^3 and 1/r^6 averaged distance below the",
528 "maximum distance ([TT]-max[tt], which will default to 0.6 in this case)",
529 "can be generated, by default averaging over equivalent hydrogens",
530 "(all triplets of hydrogens named *[123]). Additionally a list of",
531 "equivalent atoms can be supplied ([TT]-equiv[tt]), each line containing",
532 "a set of equivalent atoms specified as residue number and name and",
533 "atom name; e.g.:[PAR]",
534 "[TT]3 SER HB1 3 SER HB2[tt][PAR]",
535 "Residue and atom names must exactly match those in the structure",
536 "file, including case. Specifying non-sequential atoms is undefined."
540 int natom,i,j,teller,gi,gj;
552 atom_id *index, *noe_index;
554 real **d_r,**d,**dtot,**dtot2,**mean,**rms,**rmsc,*resnr;
555 real **dtot1_3=NULL,**dtot1_6=NULL;
556 real rmsnow,meanmax,rmsmax,rmscmax;
558 t_noe_gr *noe_gr=NULL;
562 gmx_bool bRMS, bScale, bMean, bNOE, bNMR3, bNMR6, bNMR;
564 static int nlevels=40;
565 static real scalemax=-1.0;
566 static gmx_bool bSumH=TRUE;
567 static gmx_bool bPBC=TRUE;
571 { "-nlevels", FALSE, etINT, {&nlevels},
572 "Discretize RMS in this number of levels" },
573 { "-max", FALSE, etREAL, {&scalemax},
574 "Maximum level in matrices" },
575 { "-sumh", FALSE, etBOOL, {&bSumH},
576 "Average distance over equivalent hydrogens" },
577 { "-pbc", FALSE, etBOOL, {&bPBC},
578 "Use periodic boundary conditions when computing distances" }
581 { efTRX, "-f", NULL, ffREAD },
582 { efTPS, NULL, NULL, ffREAD },
583 { efNDX, NULL, NULL, ffOPTRD },
584 { efDAT, "-equiv","equiv", ffOPTRD },
585 { efXVG, NULL, "distrmsd", ffWRITE },
586 { efXPM, "-rms", "rmsdist", ffOPTWR },
587 { efXPM, "-scl", "rmsscale", ffOPTWR },
588 { efXPM, "-mean","rmsmean", ffOPTWR },
589 { efXPM, "-nmr3","nmr3", ffOPTWR },
590 { efXPM, "-nmr6","nmr6", ffOPTWR },
591 { efDAT, "-noe", "noe", ffOPTWR },
593 #define NFILE asize(fnm)
595 CopyRight(stderr,argv[0]);
596 parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE ,
597 NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
599 bRMS = opt2bSet("-rms", NFILE,fnm);
600 bScale= opt2bSet("-scl", NFILE,fnm);
601 bMean = opt2bSet("-mean",NFILE,fnm);
602 bNOE = opt2bSet("-noe", NFILE,fnm);
603 bNMR3 = opt2bSet("-nmr3",NFILE,fnm);
604 bNMR6 = opt2bSet("-nmr6",NFILE,fnm);
605 bNMR = bNMR3 || bNMR6 || bNOE;
611 if (bNOE && scalemax < 0) {
613 fprintf(stderr,"WARNING: using -noe without -max "
614 "makes no sense, setting -max to %g\n\n",scalemax);
617 /* get topology and index */
618 read_tps_conf(ftp2fn(efTPS,NFILE,fnm),buf,&top,&ePBC,&x,NULL,box,FALSE);
624 get_index(atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&isize,&index,&grpname);
626 /* initialize arrays */
639 for(i=0; (i<isize); i++) {
642 snew(dtot2[i],isize);
644 snew(dtot1_3[i],isize);
645 snew(dtot1_6[i],isize);
655 calc_dist(isize,index,x,ePBC,box,d_r);
658 /*open output files*/
659 fp=xvgropen(ftp2fn(efXVG,NFILE,fnm),"RMS Deviation","Time (ps)","RMSD (nm)",
661 if (output_env_get_print_xvgr_codes(oenv))
662 fprintf(fp,"@ subtitle \"of distances between %s atoms\"\n",grpname);
665 natom=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
668 calc_dist_tot(isize,index,x,ePBC,box,d,dtot,dtot2,bNMR,dtot1_3,dtot1_6);
670 rmsnow=rms_diff(isize,d,d_r);
671 fprintf(fp,"%g %g\n",t,rmsnow);
672 } while (read_next_x(oenv,status,&t,natom,x,box));
673 fprintf(stderr, "\n");
677 teller = nframes_read(status);
681 calc_rms(isize,teller,dtot,dtot2,mean,&meanmax,rms,&rmsmax,rmsc,&rmscmax);
682 fprintf(stderr,"rmsmax = %g, rmscmax = %g\n",rmsmax,rmscmax);
686 calc_nmr(isize,teller,dtot1_3,dtot1_6,&max1_3,&max1_6);
689 if (scalemax > -1.0) {
698 /* make list of noe atom groups */
699 snew(noe_index, isize+1);
701 gnr=analyze_noe_equivalent(opt2fn_null("-equiv",NFILE,fnm),
702 atoms, isize, index, bSumH, noe_index, noe_gr);
703 fprintf(stdout, "Found %d non-equivalent atom-groups in %d atoms\n",
705 /* make half matrix of of noe-group distances from atom distances */
709 calc_noe(isize, noe_index, dtot1_3, dtot1_6, gnr, noe);
712 rlo.r=1.0, rlo.g=1.0, rlo.b=1.0;
713 rhi.r=0.0, rhi.g=0.0, rhi.b=0.0;
716 write_xpm(opt2FILE("-rms",NFILE,fnm,"w"),0,
717 "RMS of distance","RMS (nm)","Atom Index","Atom Index",
718 isize,isize,resnr,resnr,rms,0.0,rmsmax,rlo,rhi,&nlevels);
721 write_xpm(opt2FILE("-scl",NFILE,fnm,"w"),0,
722 "Relative RMS","RMS","Atom Index","Atom Index",
723 isize,isize,resnr,resnr,rmsc,0.0,rmscmax,rlo,rhi,&nlevels);
726 write_xpm(opt2FILE("-mean",NFILE,fnm,"w"),0,
727 "Mean Distance","Distance (nm)","Atom Index","Atom Index",
728 isize,isize,resnr,resnr,mean,0.0,meanmax,rlo,rhi,&nlevels);
731 write_xpm(opt2FILE("-nmr3",NFILE,fnm,"w"),0,"1/r^3 averaged distances",
732 "Distance (nm)","Atom Index","Atom Index",
733 isize,isize,resnr,resnr,dtot1_3,0.0,max1_3,rlo,rhi,&nlevels);
735 write_xpm(opt2FILE("-nmr6",NFILE,fnm,"w"),0,"1/r^6 averaged distances",
736 "Distance (nm)","Atom Index","Atom Index",
737 isize,isize,resnr,resnr,dtot1_6,0.0,max1_6,rlo,rhi,&nlevels);
740 write_noe(opt2FILE("-noe",NFILE,fnm,"w"), gnr, noe, noe_gr, scalemax);
742 do_view(oenv,ftp2fn(efXVG,NFILE,fnm),NULL);