4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Copyright (c) 1991-2001, University of Groningen, The Netherlands
12 * This program is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU General Public License
14 * as published by the Free Software Foundation; either version 2
15 * of the License, or (at your option) any later version.
17 * If you want to redistribute modifications, please consider that
18 * scientific software is very special. Version control is crucial -
19 * bugs must be traceable. We will be happy to consider code for
20 * inclusion in the official distribution, but derived work must not
21 * be called official GROMACS. Details are found in the README & COPYING
22 * files - if they are missing, get the official version at www.gromacs.org.
24 * To help us fund GROMACS development, we humbly ask that you cite
25 * the papers on the package - you can find them in the top README file.
27 * For more info, check our website at http://www.gromacs.org
30 * GROtesk MACabre and Sinister
32 static char *SRCID_g_rmsdist_c
= "$Id$";
53 static void calc_dist(int nind
,atom_id index
[],rvec x
[],real
**d
)
59 for(i
=0; (i
<nind
-1); i
++) {
61 for(j
=i
+1; (j
<nind
); j
++) {
62 pbc_dx(xi
,x
[index
[j
]],dx
);
68 static void calc_dist_tot(int nind
,atom_id index
[], rvec x
[],
69 real
**d
, real
**dtot
, real
**dtot2
,
70 bool bNMR
, real
**dtot1_3
, real
**dtot1_6
)
74 real temp
, temp2
, temp1_3
;
77 for(i
=0; (i
<nind
-1); i
++) {
79 for(j
=i
+1; (j
<nind
); j
++) {
80 pbc_dx(xi
,x
[index
[j
]],dx
);
81 temp2
=dx
[XX
]*dx
[XX
]+dx
[YY
]*dx
[YY
]+dx
[ZZ
]*dx
[ZZ
];
87 temp1_3
= 1.0/(temp
*temp2
);
88 dtot1_3
[i
][j
] += temp1_3
;
89 dtot1_6
[i
][j
] += temp1_3
*temp1_3
;
95 static void calc_nmr(int nind
, int nframes
, real
**dtot1_3
, real
**dtot1_6
,
96 real
*max1_3
, real
*max1_6
)
101 for(i
=0; (i
<nind
-1); i
++) {
102 for(j
=i
+1; (j
<nind
); j
++) {
103 temp1_3
= pow(dtot1_3
[i
][j
]/nframes
,-1.0/3.0);
104 temp1_6
= pow(dtot1_6
[i
][j
]/nframes
,-1.0/6.0);
105 if (temp1_3
> *max1_3
) *max1_3
= temp1_3
;
106 if (temp1_6
> *max1_6
) *max1_6
= temp1_6
;
107 dtot1_3
[i
][j
] = temp1_3
;
108 dtot1_6
[i
][j
] = temp1_6
;
109 dtot1_3
[j
][i
] = temp1_3
;
110 dtot1_6
[j
][i
] = temp1_6
;
115 static char Hnum
[] = "123";
140 static int read_equiv(char *eq_fn
, t_equiv
***equivptr
)
143 char line
[STRLEN
],resname
[10],atomname
[10],*lp
;
144 int neq
, na
, n
, resnr
;
147 fp
= ffopen(eq_fn
,"r");
150 while(get_a_line(fp
,line
,STRLEN
)) {
152 /* this is not efficient, but I'm lazy */
156 if (sscanf(lp
,"%s %n",atomname
,&n
)==1) {
159 equiv
[neq
][0].nname
=strdup(atomname
);
160 while (sscanf(lp
,"%d %s %s %n",&resnr
,resname
,atomname
,&n
)==3) {
161 /* this is not efficient, but I'm lazy (again) */
162 srenew(equiv
[neq
], na
+1);
163 equiv
[neq
][na
].rnr
=resnr
-1;
164 equiv
[neq
][na
].rname
=strdup(resname
);
165 equiv
[neq
][na
].aname
=strdup(atomname
);
166 if (na
>0) equiv
[neq
][na
].nname
=NULL
;
171 /* make empty element as flag for end of array */
172 srenew(equiv
[neq
], na
+1);
173 equiv
[neq
][na
].rnr
=NOTSET
;
174 equiv
[neq
][na
].rname
=NULL
;
175 equiv
[neq
][na
].aname
=NULL
;
187 static void dump_equiv(FILE *out
, int neq
, t_equiv
**equiv
)
191 fprintf(out
,"Dumping equivalent list\n");
192 for (i
=0; i
<neq
; i
++) {
193 fprintf(out
,"%s",equiv
[i
][0].nname
);
194 for(j
=0; equiv
[i
][j
].rnr
!=NOTSET
; j
++)
195 fprintf(out
," %d %s %s",
196 equiv
[i
][j
].rnr
,equiv
[i
][j
].rname
,equiv
[i
][j
].aname
);
201 static bool is_equiv(int neq
, t_equiv
**equiv
, char **nname
,
202 int rnr1
, char *rname1
, char *aname1
,
203 int rnr2
, char *rname2
, char *aname2
)
209 /* we can terminate each loop when bFound is true! */
210 for (i
=0; i
<neq
&& !bFound
; i
++) {
211 /* find first atom */
212 for(j
=0; equiv
[i
][j
].rnr
!=NOTSET
&& !bFound
; j
++)
213 bFound
= ( equiv
[i
][j
].rnr
==rnr1
&&
214 strcmp(equiv
[i
][j
].rname
, rname1
)==0 &&
215 strcmp(equiv
[i
][j
].aname
, aname1
)==0 );
217 /* find second atom */
219 for(j
=0; equiv
[i
][j
].rnr
!=NOTSET
&& !bFound
; j
++)
220 bFound
= ( equiv
[i
][j
].rnr
==rnr2
&&
221 strcmp(equiv
[i
][j
].rname
, rname2
)==0 &&
222 strcmp(equiv
[i
][j
].aname
, aname2
)==0 );
226 *nname
= strdup(equiv
[i
-1][0].nname
);
231 static int analyze_noe_equivalent(char *eq_fn
,
232 t_atoms
*atoms
, int isize
, atom_id
*index
,
234 atom_id
*noe_index
, t_noe_gr
*noe_gr
)
237 int i
, j
, anmil
, anmjl
, rnri
, rnrj
, gi
, groupnr
, neq
;
238 char *anmi
, *anmj
, **nnm
;
245 neq
=read_equiv(eq_fn
,&equiv
);
246 if (debug
) dump_equiv(debug
,neq
,equiv
);
253 for(i
=0; i
<isize
; i
++) {
254 if (equiv
&& i
<isize
-1) {
255 /* check explicit list of equivalent atoms */
258 rnri
=atoms
->atom
[index
[i
]].resnr
;
259 rnrj
=atoms
->atom
[index
[j
]].resnr
;
261 is_equiv(neq
, equiv
, &nnm
[i
],
262 rnri
, *atoms
->resname
[rnri
], *atoms
->atomname
[index
[i
]],
263 rnrj
, *atoms
->resname
[rnrj
], *atoms
->atomname
[index
[j
]]);
265 nnm
[j
]=strdup(nnm
[i
]);
267 /* set index for matching atom */
268 noe_index
[j
]=groupnr
;
269 /* skip matching atom */
272 } while ( bEquiv
&& i
<isize
-1 );
276 /* look for triplets of consecutive atoms with name XX?,
277 X are any number of letters or digits and ? goes from 1 to 3
278 This is supposed to cover all CH3 groups and the like */
279 anmi
= *atoms
->atomname
[index
[i
]];
280 anmil
= strlen(anmi
);
281 bMatch
= i
<isize
-3 && anmi
[anmil
-1]=='1';
284 anmj
= *atoms
->atomname
[index
[i
+j
]];
285 anmjl
= strlen(anmj
);
286 bMatch
= bMatch
&& ( anmil
==anmjl
&& anmj
[anmjl
-1]==Hnum
[j
] &&
287 strncmp(anmi
, anmj
, anmil
-1)==0 );
289 /* set index for this atom */
290 noe_index
[i
]=groupnr
;
292 /* set index for next two matching atoms */
294 noe_index
[i
+j
]=groupnr
;
295 /* skip matching atoms */
302 /* make index without looking for equivalent atoms */
303 for(i
=0; i
<isize
; i
++)
307 noe_index
[isize
]=groupnr
;
311 for(i
=0; i
<isize
; i
++) {
312 rnri
=atoms
->atom
[index
[i
]].resnr
;
313 fprintf(debug
,"%s %s %d -> %s\n",*atoms
->atomname
[index
[i
]],
314 *atoms
->resname
[rnri
],rnri
,nnm
[i
]?nnm
[i
]:"");
317 for(i
=0; i
<isize
; i
++) {
319 if (!noe_gr
[gi
].aname
) {
321 noe_gr
[gi
].anr
=index
[i
];
323 noe_gr
[gi
].aname
=strdup(nnm
[i
]);
325 noe_gr
[gi
].aname
=strdup(*atoms
->atomname
[index
[i
]]);
326 if ( noe_index
[i
]==noe_index
[i
+1] )
327 noe_gr
[gi
].aname
[strlen(noe_gr
[gi
].aname
)-1]='*';
329 noe_gr
[gi
].rnr
=atoms
->atom
[index
[i
]].resnr
;
330 noe_gr
[gi
].rname
=strdup(*atoms
->resname
[noe_gr
[gi
].rnr
]);
331 /* dump group definitions */
332 if (debug
) fprintf(debug
,"%d %d %d %d %s %s %d\n",i
,gi
,
333 noe_gr
[gi
].ianr
,noe_gr
[gi
].anr
,noe_gr
[gi
].aname
,
334 noe_gr
[gi
].rname
,noe_gr
[gi
].rnr
);
337 for(i
=0; i
<isize
; i
++)
344 /* #define NSCALE 3 */
345 /* static char *noe_scale[NSCALE+1] = { */
346 /* "strong", "medium", "weak", "none" */
350 static char *noe2scale(real r3
, real r6
, real rmax
)
353 static char buf
[NSCALE
+1];
355 /* r goes from 0 to rmax
356 NSCALE*r/rmax goes from 0 to NSCALE
357 NSCALE - NSCALE*r/rmax goes from NSCALE to 0 */
358 s3
= NSCALE
- min(NSCALE
, (int)(NSCALE
*r3
/rmax
));
359 s6
= NSCALE
- min(NSCALE
, (int)(NSCALE
*r6
/rmax
));
370 static void calc_noe(int isize
, atom_id
*noe_index
,
371 real
**dtot1_3
, real
**dtot1_6
, int gnr
, t_noe
**noe
)
375 /* make half matrix of noe-group distances from atom distances */
376 for(i
=0; i
<isize
; i
++) {
378 for(j
=i
; j
<isize
; j
++) {
381 noe
[gi
][gj
].i_3
+=pow(dtot1_3
[i
][j
],-3);
382 noe
[gi
][gj
].i_6
+=pow(dtot1_6
[i
][j
],-6);
388 for(j
=i
+1; j
<gnr
; j
++) {
389 noe
[i
][j
].r_3
= pow(noe
[i
][j
].i_3
/noe
[i
][j
].nr
,-1.0/3.0);
390 noe
[i
][j
].r_6
= pow(noe
[i
][j
].i_6
/noe
[i
][j
].nr
,-1.0/6.0);
391 noe
[j
][i
] = noe
[i
][j
];
395 static void write_noe(FILE *fp
,int gnr
,t_noe
**noe
,t_noe_gr
*noe_gr
,real rmax
)
398 real r3
, r6
, min3
, min6
;
399 char buf
[10],b3
[10],b6
[10];
404 ";%4s %3s %4s %4s%3s %4s %4s %4s %4s%3s %5s %5s %8s %2s %2s %s\n",
405 "ianr","anr","anm","rnm","rnr","ianr","anr","anm","rnm","rnr",
406 "1/r^3","1/r^6","intnsty","Dr","Da","scale");
407 for(i
=0; i
<gnr
; i
++) {
409 for(j
=i
+1; j
<gnr
; j
++) {
415 if ( r3
< rmax
|| r6
< rmax
) {
416 if (grj
.rnr
== gri
.rnr
)
417 sprintf(buf
,"%2d", grj
.anr
-gri
.anr
);
421 sprintf(b3
,"%-5.3f",r3
);
425 sprintf(b6
,"%-5.3f",r6
);
429 "%4d %4d %4s %4s%3d %4d %4d %4s %4s%3d %5s %5s %8d %2d %2s %s\n",
430 gri
.ianr
+1, gri
.anr
+1, gri
.aname
, gri
.rname
, gri
.rnr
+1,
431 grj
.ianr
+1, grj
.anr
+1, grj
.aname
, grj
.rname
, grj
.rnr
+1,
432 b3
, b6
, (int)(noe
[i
][j
].i_6
+0.5), grj
.rnr
-gri
.rnr
, buf
,
433 noe2scale(r3
, r6
, rmax
));
437 #define MINI ((i==3)?min3:min6)
440 fprintf(stdout
,"NOTE: no 1/r^%d averaged distances found below %g, "
441 "smallest was %g\n", i
, rmax
, MINI
);
443 fprintf(stdout
,"Smallest 1/r^%d averaged distance was %g\n", i
, MINI
);
447 static void calc_rms(int nind
, int nframes
,
448 real
**dtot
, real
**dtot2
,
449 real
**rmsmat
, real
*rmsmax
,
450 real
**rmscmat
, real
*rmscmax
,
451 real
**meanmat
, real
*meanmax
)
454 real mean
, mean2
, rms
, rmsc
;
455 /* N.B. dtot and dtot2 contain the total distance and the total squared
456 * distance respectively, BUT they return RMS and the scaled RMS resp.
463 for(i
=0; (i
<nind
-1); i
++) {
464 for(j
=i
+1; (j
<nind
); j
++) {
465 mean
=dtot
[i
][j
]/nframes
;
466 mean2
=dtot2
[i
][j
]/nframes
;
467 rms
=sqrt(max(0,mean2
-mean
*mean
));
469 if (mean
> *meanmax
) *meanmax
=mean
;
470 if (rms
> *rmsmax
) *rmsmax
=rms
;
471 if (rmsc
> *rmscmax
) *rmscmax
=rmsc
;
472 meanmat
[i
][j
]=meanmat
[j
][i
]=mean
;
473 rmsmat
[i
][j
] =rmsmat
[j
][i
] =rms
;
474 rmscmat
[i
][j
]=rmscmat
[j
][i
]=rmsc
;
479 real
rms_diff(int natom
,real
**d
,real
**d_r
)
485 for(i
=0; (i
<natom
-1); i
++)
486 for(j
=i
+1; (j
<natom
); j
++) {
490 r2
/=(natom
*(natom
-1))/2;
495 int main (int argc
,char *argv
[])
497 static char *desc
[] = {
498 "g_rmsdist computes the root mean square deviation of atom distances,",
499 "which has the advantage that no fit is needed like in standard RMS",
500 "deviation as computed by g_rms.",
501 "The reference structure is taken from the structure file.",
502 "The rmsd at time t is calculated as the rms",
503 "of the differences in distance between atom-pairs in the reference",
504 "structure and the structure at time t.[PAR]",
505 "g_rmsdist can also produce matrices of the rms distances, rms distances",
506 "scaled with the mean distance and the mean distances and matrices with",
507 "NMR averaged distances (1/r^3 and 1/r^6 averaging). Finally, lists",
508 "of atom pairs with 1/r^3 and 1/r^6 averaged distance below the",
509 "maximum distance ([TT]-max[tt], which will default to 0.6 in this case)",
510 "can be generated, by default averaging over equivalent hydrogens",
511 "(all triplets of hydrogens named *[123]). Additionally a list of",
512 "equivalent atoms can be supplied ([TT]-equiv[tt]), each line containing",
513 "a set of equivalent atoms specified as residue number and name and",
514 "atom name; e.g.:[PAR]",
515 "[TT]3 SER HB1 3 SER HB2[tt][PAR]",
516 "Residue and atom names must exactly match those in the structure",
517 "file, including case. Specifying non-sequential atoms is undefined."
521 int natom
,i
,j
,teller
,gi
,gj
;
530 int status
,isize
,gnr
=0;
531 atom_id
*index
, *noe_index
;
533 real
**d_r
,**d
,**dtot
,**dtot2
,**mean
,**rms
,**rmsc
,*resnr
;
534 real
**dtot1_3
=NULL
,**dtot1_6
=NULL
;
535 real rmsnow
,meanmax
,rmsmax
,rmscmax
;
537 t_noe_gr
*noe_gr
=NULL
;
541 bool bRMS
, bScale
, bMean
, bNOE
, bNMR3
, bNMR6
, bNMR
;
543 static int nlevels
=40;
544 static real scalemax
=-1.0;
545 static bool bSumH
=TRUE
;
547 { "-nlevels", FALSE
, etINT
, {&nlevels
},
548 "Discretize rms in # levels" },
549 { "-max", FALSE
, etREAL
, {&scalemax
},
550 "Maximum level in matrices" },
551 { "-sumh", FALSE
, etBOOL
, {&bSumH
},
552 "average distance over equivalent hydrogens" }
555 { efTRX
, "-f", NULL
, ffREAD
},
556 { efTPS
, NULL
, NULL
, ffREAD
},
557 { efNDX
, NULL
, NULL
, ffOPTRD
},
558 { efDAT
, "-equiv","equiv", ffOPTRD
},
559 { efXVG
, NULL
, "distrmsd", ffWRITE
},
560 { efXPM
, "-rms", "rmsdist", ffOPTWR
},
561 { efXPM
, "-scl", "rmsscale", ffOPTWR
},
562 { efXPM
, "-mean","rmsmean", ffOPTWR
},
563 { efXPM
, "-nmr3","nmr3", ffOPTWR
},
564 { efXPM
, "-nmr6","nmr6", ffOPTWR
},
565 { efDAT
, "-noe", "noe", ffOPTWR
},
567 #define NFILE asize(fnm)
569 CopyRight(stderr
,argv
[0]);
570 parse_common_args(&argc
,argv
,PCA_CAN_VIEW
| PCA_CAN_TIME
| PCA_BE_NICE
,
571 NFILE
,fnm
,asize(pa
),pa
,asize(desc
),desc
,0,NULL
);
573 bRMS
= opt2bSet("-rms", NFILE
,fnm
);
574 bScale
= opt2bSet("-scl", NFILE
,fnm
);
575 bMean
= opt2bSet("-mean",NFILE
,fnm
);
576 bNOE
= opt2bSet("-noe", NFILE
,fnm
);
577 bNMR3
= opt2bSet("-nmr3",NFILE
,fnm
);
578 bNMR6
= opt2bSet("-nmr6",NFILE
,fnm
);
579 bNMR
= bNMR3
|| bNMR6
|| bNOE
;
582 if (bNOE
&& scalemax
< 0) {
584 fprintf(stderr
,"WARNING: using -noe without -max "
585 "makes no sense, setting -max to %g\n\n",scalemax
);
588 /* get topology and index */
589 read_tps_conf(ftp2fn(efTPS
,NFILE
,fnm
),buf
,&top
,&x
,NULL
,box
,FALSE
);
592 get_index(atoms
,ftp2fn_null(efNDX
,NFILE
,fnm
),1,&isize
,&index
,&grpname
);
594 /* initialize arrays */
607 for(i
=0; (i
<isize
); i
++) {
610 snew(dtot2
[i
],isize
);
612 snew(dtot1_3
[i
],isize
);
613 snew(dtot1_6
[i
],isize
);
624 calc_dist(isize
,index
,x
,d_r
);
627 /*open output files*/
628 fp
=xvgropen(ftp2fn(efXVG
,NFILE
,fnm
),"RMS Deviation","Time (ps)","RMSD (nm)");
629 fprintf(fp
,"@ subtitle \"of distances between %s atoms\"\n",grpname
);
632 natom
=read_first_x(&status
,ftp2fn(efTRX
,NFILE
,fnm
),&t
,&x
,box
);
635 calc_dist_tot(isize
,index
,x
,d
,dtot
,dtot2
,bNMR
,dtot1_3
,dtot1_6
);
637 rmsnow
=rms_diff(isize
,d
,d_r
);
638 fprintf(fp
,"%g %g\n",t
,rmsnow
);
639 } while (read_next_x(status
,&t
,natom
,x
,box
));
640 fprintf(stderr
, "\n");
645 teller
= nframes_read();
646 calc_rms(isize
,teller
,dtot
,dtot2
,mean
,&meanmax
,rms
,&rmsmax
,rmsc
,&rmscmax
);
647 fprintf(stderr
,"rmsmax = %g, rmscmax = %g\n",rmsmax
,rmscmax
);
652 calc_nmr(isize
,teller
,dtot1_3
,dtot1_6
,&max1_3
,&max1_6
);
655 if (scalemax
> -1.0) {
664 /* make list of noe atom groups */
665 snew(noe_index
, isize
+1);
667 gnr
=analyze_noe_equivalent(opt2fn_null("-equiv",NFILE
,fnm
),
668 atoms
, isize
, index
, bSumH
, noe_index
, noe_gr
);
669 fprintf(stdout
, "Found %d non-equivalent atom-groups in %d atoms\n",
671 /* make half matrix of of noe-group distances from atom distances */
675 calc_noe(isize
, noe_index
, dtot1_3
, dtot1_6
, gnr
, noe
);
678 rlo
.r
=1.0, rlo
.g
=1.0, rlo
.b
=1.0;
679 rhi
.r
=0.0, rhi
.g
=0.0, rhi
.b
=0.0;
682 write_xpm(opt2FILE("-rms",NFILE
,fnm
,"w"),
683 "RMS of distance","RMS (nm)","Atom Index","Atom Index",
684 isize
,isize
,resnr
,resnr
,rms
,0.0,rmsmax
,rlo
,rhi
,&nlevels
);
687 write_xpm(opt2FILE("-scl",NFILE
,fnm
,"w"),
688 "Relative RMS","RMS","Atom Index","Atom Index",
689 isize
,isize
,resnr
,resnr
,rmsc
,0.0,rmscmax
,rlo
,rhi
,&nlevels
);
692 write_xpm(opt2FILE("-mean",NFILE
,fnm
,"w"),
693 "Mean Distance","Distance (nm)","Atom Index","Atom Index",
694 isize
,isize
,resnr
,resnr
,mean
,0.0,meanmax
,rlo
,rhi
,&nlevels
);
697 write_xpm(opt2FILE("-nmr3",NFILE
,fnm
,"w"),"1/r^3 averaged distances",
698 "Distance (nm)","Atom Index","Atom Index",
699 isize
,isize
,resnr
,resnr
,dtot1_3
,0.0,max1_3
,rlo
,rhi
,&nlevels
);
701 write_xpm(opt2FILE("-nmr6",NFILE
,fnm
,"w"),"1/r^6 averaged distances",
702 "Distance (nm)","Atom Index","Atom Index",
703 isize
,isize
,resnr
,resnr
,dtot1_6
,0.0,max1_6
,rlo
,rhi
,&nlevels
);
706 write_noe(opt2FILE("-noe",NFILE
,fnm
,"w"), gnr
, noe
, noe_gr
, scalemax
);
708 do_view(ftp2fn(efXVG
,NFILE
,fnm
),NULL
);