?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Upped the version to 3.2.0
[gromacs.git] / src / tools / g_rmsdist.c
blob041f57af7296b8c3c206a230f7d058d684575d8b
1 /*
2 * $Id$
3 *
4 * This source code is part of
5 *
6 * G R O M A C S
7 *
8 * GROningen MAchine for Chemical Simulations
9 *
10 * VERSION 3.2.0
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
15
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
20 *
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
27 *
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
30 *
31 * For more info, check our website at http://www.gromacs.org
32 *
33 * And Hey:
34 * Green Red Orange Magenta Azure Cyan Skyblue
35 */
36 #include <math.h>
37 #include <ctype.h>
38 #include "macros.h"
39 #include "smalloc.h"
40 #include "typedefs.h"
41 #include "names.h"
42 #include "copyrite.h"
43 #include "statutil.h"
44 #include "tpxio.h"
45 #include "string2.h"
46 #include "strdb.h"
47 #include "vec.h"
48 #include "macros.h"
49 #include "index.h"
50 #include "pbc.h"
51 #include "xvgr.h"
52 #include "futil.h"
53 #include "matio.h"
54 #include "assert.h"
55
56 static void calc_dist(int nind,atom_id index[],rvec x[],real **d)
57 {
58 int i,j;
59 real *xi;
60 rvec dx;
61
62 for(i=0; (i<nind-1); i++) {
63 xi=x[index[i]];
64 for(j=i+1; (j<nind); j++) {
65 pbc_dx(xi,x[index[j]],dx);
66 d[i][j]=norm(dx);
67 }
68 }
69 }
70
71 static void calc_dist_tot(int nind,atom_id index[], rvec x[],
72 real **d, real **dtot, real **dtot2,
73 bool bNMR, real **dtot1_3, real **dtot1_6)
74 {
75 int i,j;
76 real *xi;
77 real temp, temp2, temp1_3;
78 rvec dx;
79
80 for(i=0; (i<nind-1); i++) {
81 xi=x[index[i]];
82 for(j=i+1; (j<nind); j++) {
83 pbc_dx(xi,x[index[j]],dx);
84 temp2=dx[XX]*dx[XX]+dx[YY]*dx[YY]+dx[ZZ]*dx[ZZ];
85 temp =sqrt(temp2);
86 d[i][j]=temp;
87 dtot[i][j]+=temp;
88 dtot2[i][j]+=temp2;
89 if (bNMR) {
90 temp1_3 = 1.0/(temp*temp2);
91 dtot1_3[i][j] += temp1_3;
92 dtot1_6[i][j] += temp1_3*temp1_3;
93 }
94 }
95 }
96 }
97
98 static void calc_nmr(int nind, int nframes, real **dtot1_3, real **dtot1_6,
99 real *max1_3, real *max1_6)
100 {
101 int i,j;
102 real temp1_3,temp1_6;
103
104 for(i=0; (i<nind-1); i++) {
105 for(j=i+1; (j<nind); j++) {
106 temp1_3 = pow(dtot1_3[i][j]/nframes,-1.0/3.0);
107 temp1_6 = pow(dtot1_6[i][j]/nframes,-1.0/6.0);
108 if (temp1_3 > *max1_3) *max1_3 = temp1_3;
109 if (temp1_6 > *max1_6) *max1_6 = temp1_6;
110 dtot1_3[i][j] = temp1_3;
111 dtot1_6[i][j] = temp1_6;
112 dtot1_3[j][i] = temp1_3;
113 dtot1_6[j][i] = temp1_6;
114 }
115 }
116 }
117
118 static char Hnum[] = "123";
119
120 typedef struct {
121 int nr;
122 real r_3;
123 real r_6;
124 real i_3;
125 real i_6;
126 } t_noe;
127
128 typedef struct {
129 int anr;
130 int ianr;
131 int rnr;
132 char *aname;
133 char *rname;
134 } t_noe_gr;
135
136 typedef struct {
137 int rnr;
138 char *nname;
139 char *rname;
140 char *aname;
141 } t_equiv;
142
143 static int read_equiv(char *eq_fn, t_equiv ***equivptr)
144 {
145 FILE *fp;
146 char line[STRLEN],resname[10],atomname[10],*lp;
147 int neq, na, n, resnr;
148 t_equiv **equiv;
149
150 fp = ffopen(eq_fn,"r");
151 neq=0;
152 equiv=NULL;
153 while(get_a_line(fp,line,STRLEN)) {
154 lp=line;
155 /* this is not efficient, but I'm lazy */
156 srenew(equiv,neq+1);
157 equiv[neq]=NULL;
158 na=0;
159 if (sscanf(lp,"%s %n",atomname,&n)==1) {
160 lp+=n;
161 snew(equiv[neq], 1);
162 equiv[neq][0].nname=strdup(atomname);
163 while (sscanf(lp,"%d %s %s %n",&resnr,resname,atomname,&n)==3) {
164 /* this is not efficient, but I'm lazy (again) */
165 srenew(equiv[neq], na+1);
166 equiv[neq][na].rnr=resnr-1;
167 equiv[neq][na].rname=strdup(resname);
168 equiv[neq][na].aname=strdup(atomname);
169 if (na>0) equiv[neq][na].nname=NULL;
170 na++;
171 lp+=n;
172 }
173 }
174 /* make empty element as flag for end of array */
175 srenew(equiv[neq], na+1);
176 equiv[neq][na].rnr=NOTSET;
177 equiv[neq][na].rname=NULL;
178 equiv[neq][na].aname=NULL;
179
180 /* next */
181 neq++;
182 }
183 ffclose(fp);
184
185 *equivptr=equiv;
186
187 return neq;
188 }
189
190 static void dump_equiv(FILE *out, int neq, t_equiv **equiv)
191 {
192 int i,j;
193
194 fprintf(out,"Dumping equivalent list\n");
195 for (i=0; i<neq; i++) {
196 fprintf(out,"%s",equiv[i][0].nname);
197 for(j=0; equiv[i][j].rnr!=NOTSET; j++)
198 fprintf(out," %d %s %s",
199 equiv[i][j].rnr,equiv[i][j].rname,equiv[i][j].aname);
200 fprintf(out,"\n");
201 }
202 }
203
204 static bool is_equiv(int neq, t_equiv **equiv, char **nname,
205 int rnr1, char *rname1, char *aname1,
206 int rnr2, char *rname2, char *aname2)
207 {
208 int i,j;
209 bool bFound;
210
211 bFound=FALSE;
212 /* we can terminate each loop when bFound is true! */
213 for (i=0; i<neq && !bFound; i++) {
214 /* find first atom */
215 for(j=0; equiv[i][j].rnr!=NOTSET && !bFound; j++)
216 bFound = ( equiv[i][j].rnr==rnr1 &&
217 strcmp(equiv[i][j].rname, rname1)==0 &&
218 strcmp(equiv[i][j].aname, aname1)==0 );
219 if (bFound) {
220 /* find second atom */
221 bFound=FALSE;
222 for(j=0; equiv[i][j].rnr!=NOTSET && !bFound; j++)
223 bFound = ( equiv[i][j].rnr==rnr2 &&
224 strcmp(equiv[i][j].rname, rname2)==0 &&
225 strcmp(equiv[i][j].aname, aname2)==0 );
226 }
227 }
228 if (bFound)
229 *nname = strdup(equiv[i-1][0].nname);
230
231 return bFound;
232 }
233
234 static int analyze_noe_equivalent(char *eq_fn,
235 t_atoms *atoms, int isize, atom_id *index,
236 bool bSumH,
237 atom_id *noe_index, t_noe_gr *noe_gr)
238 {
239 FILE *fp;
240 int i, j, anmil, anmjl, rnri, rnrj, gi, groupnr, neq;
241 char *anmi, *anmj, **nnm;
242 bool bMatch,bEquiv;
243 t_equiv **equiv;
244
245 snew(nnm,isize);
246 if (bSumH) {
247 if (eq_fn) {
248 neq=read_equiv(eq_fn,&equiv);
249 if (debug) dump_equiv(debug,neq,equiv);
250 } else {
251 neq=0;
252 equiv=NULL;
253 }
254
255 groupnr=0;
256 for(i=0; i<isize; i++) {
257 if (equiv && i<isize-1) {
258 /* check explicit list of equivalent atoms */
259 do {
260 j=i+1;
261 rnri=atoms->atom[index[i]].resnr;
262 rnrj=atoms->atom[index[j]].resnr;
263 bEquiv =
264 is_equiv(neq, equiv, &nnm[i],
265 rnri, *atoms->resname[rnri], *atoms->atomname[index[i]],
266 rnrj, *atoms->resname[rnrj], *atoms->atomname[index[j]]);
267 if(nnm[i] && bEquiv)
268 nnm[j]=strdup(nnm[i]);
269 if (bEquiv) {
270 /* set index for matching atom */
271 noe_index[j]=groupnr;
272 /* skip matching atom */
273 i=j;
274 }
275 } while ( bEquiv && i<isize-1 );
276 } else
277 bEquiv=FALSE;
278 if (!bEquiv) {
279 /* look for triplets of consecutive atoms with name XX?,
280 X are any number of letters or digits and ? goes from 1 to 3
281 This is supposed to cover all CH3 groups and the like */
282 anmi = *atoms->atomname[index[i]];
283 anmil = strlen(anmi);
284 bMatch = i<isize-3 && anmi[anmil-1]=='1';
285 if (bMatch)
286 for(j=1; j<3; j++) {
287 anmj = *atoms->atomname[index[i+j]];
288 anmjl = strlen(anmj);
289 bMatch = bMatch && ( anmil==anmjl && anmj[anmjl-1]==Hnum[j] &&
290 strncmp(anmi, anmj, anmil-1)==0 );
291 }
292 /* set index for this atom */
293 noe_index[i]=groupnr;
294 if (bMatch) {
295 /* set index for next two matching atoms */
296 for(j=1; j<3; j++)
297 noe_index[i+j]=groupnr;
298 /* skip matching atoms */
299 i+=2;
300 }
301 }
302 groupnr++;
303 }
304 } else {
305 /* make index without looking for equivalent atoms */
306 for(i=0; i<isize; i++)
307 noe_index[i]=i;
308 groupnr=isize;
309 }
310 noe_index[isize]=groupnr;
311
312 if (debug)
313 /* dump new names */
314 for(i=0; i<isize; i++) {
315 rnri=atoms->atom[index[i]].resnr;
316 fprintf(debug,"%s %s %d -> %s\n",*atoms->atomname[index[i]],
317 *atoms->resname[rnri],rnri,nnm[i]?nnm[i]:"");
318 }
319
320 for(i=0; i<isize; i++) {
321 gi=noe_index[i];
322 if (!noe_gr[gi].aname) {
323 noe_gr[gi].ianr=i;
324 noe_gr[gi].anr=index[i];
325 if (nnm[i])
326 noe_gr[gi].aname=strdup(nnm[i]);
327 else {
328 noe_gr[gi].aname=strdup(*atoms->atomname[index[i]]);
329 if ( noe_index[i]==noe_index[i+1] )
330 noe_gr[gi].aname[strlen(noe_gr[gi].aname)-1]='*';
331 }
332 noe_gr[gi].rnr=atoms->atom[index[i]].resnr;
333 noe_gr[gi].rname=strdup(*atoms->resname[noe_gr[gi].rnr]);
334 /* dump group definitions */
335 if (debug) fprintf(debug,"%d %d %d %d %s %s %d\n",i,gi,
336 noe_gr[gi].ianr,noe_gr[gi].anr,noe_gr[gi].aname,
337 noe_gr[gi].rname,noe_gr[gi].rnr);
338 }
339 }
340 for(i=0; i<isize; i++)
341 sfree(nnm[i]);
342 sfree(nnm);
343
344 return groupnr;
345 }
346
347 /* #define NSCALE 3 */
348 /* static char *noe_scale[NSCALE+1] = { */
349 /* "strong", "medium", "weak", "none" */
350 /* }; */
351 #define NSCALE 6
352
353 static char *noe2scale(real r3, real r6, real rmax)
354 {
355 int i,s3, s6;
356 static char buf[NSCALE+1];
357
358 /* r goes from 0 to rmax
359 NSCALE*r/rmax goes from 0 to NSCALE
360 NSCALE - NSCALE*r/rmax goes from NSCALE to 0 */
361 s3 = NSCALE - min(NSCALE, (int)(NSCALE*r3/rmax));
362 s6 = NSCALE - min(NSCALE, (int)(NSCALE*r6/rmax));
363
364 for(i=0; i<s3; i++)
365 buf[i]='=';
366 for( ; i<s6; i++)
367 buf[i]='-';
368 buf[i]='\0';
369
370 return buf;
371 }
372
373 static void calc_noe(int isize, atom_id *noe_index,
374 real **dtot1_3, real **dtot1_6, int gnr, t_noe **noe)
375 {
376 int i, j, gi, gj;
377
378 /* make half matrix of noe-group distances from atom distances */
379 for(i=0; i<isize; i++) {
380 gi=noe_index[i];
381 for(j=i; j<isize; j++) {
382 gj=noe_index[j];
383 noe[gi][gj].nr++;
384 noe[gi][gj].i_3+=pow(dtot1_3[i][j],-3);
385 noe[gi][gj].i_6+=pow(dtot1_6[i][j],-6);
386 }
387 }
388
389 /* make averages */
390 for(i=0; i<gnr; i++)
391 for(j=i+1; j<gnr; j++) {
392 noe[i][j].r_3 = pow(noe[i][j].i_3/noe[i][j].nr,-1.0/3.0);
393 noe[i][j].r_6 = pow(noe[i][j].i_6/noe[i][j].nr,-1.0/6.0);
394 noe[j][i] = noe[i][j];
395 }
396 }
397
398 static void write_noe(FILE *fp,int gnr,t_noe **noe,t_noe_gr *noe_gr,real rmax)
399 {
400 int i,j;
401 real r3, r6, min3, min6;
402 char buf[10],b3[10],b6[10];
403 t_noe_gr gri, grj;
404
405 min3 = min6 = 1e6;
406 fprintf(fp,
407 ";%4s %3s %4s %4s%3s %4s %4s %4s %4s%3s %5s %5s %8s %2s %2s %s\n",
408 "ianr","anr","anm","rnm","rnr","ianr","anr","anm","rnm","rnr",
409 "1/r^3","1/r^6","intnsty","Dr","Da","scale");
410 for(i=0; i<gnr; i++) {
411 gri=noe_gr[i];
412 for(j=i+1; j<gnr; j++) {
413 grj=noe_gr[j];
414 r3 = noe[i][j].r_3;
415 r6 = noe[i][j].r_6;
416 min3 = min(r3,min3);
417 min6 = min(r6,min6);
418 if ( r3 < rmax || r6 < rmax ) {
419 if (grj.rnr == gri.rnr)
420 sprintf(buf,"%2d", grj.anr-gri.anr);
421 else
422 buf[0]='\0';
423 if ( r3 < rmax )
424 sprintf(b3,"%-5.3f",r3);
425 else
426 strcpy(b3,"-");
427 if ( r6 < rmax )
428 sprintf(b6,"%-5.3f",r6);
429 else
430 strcpy(b6,"-");
431 fprintf(fp,
432 "%4d %4d %4s %4s%3d %4d %4d %4s %4s%3d %5s %5s %8d %2d %2s %s\n",
433 gri.ianr+1, gri.anr+1, gri.aname, gri.rname, gri.rnr+1,
434 grj.ianr+1, grj.anr+1, grj.aname, grj.rname, grj.rnr+1,
435 b3, b6, (int)(noe[i][j].i_6+0.5), grj.rnr-gri.rnr, buf,
436 noe2scale(r3, r6, rmax));
437 }
438 }
439 }
440 #define MINI ((i==3)?min3:min6)
441 for(i=3; i<=6; i+=3)
442 if ( MINI > rmax )
443 fprintf(stdout,"NOTE: no 1/r^%d averaged distances found below %g, "
444 "smallest was %g\n", i, rmax, MINI );
445 else
446 fprintf(stdout,"Smallest 1/r^%d averaged distance was %g\n", i, MINI );
447 #undef MINI
448 }
449
450 static void calc_rms(int nind, int nframes,
451 real **dtot, real **dtot2,
452 real **rmsmat, real *rmsmax,
453 real **rmscmat, real *rmscmax,
454 real **meanmat, real *meanmax)
455 {
456 int i,j;
457 real mean, mean2, rms, rmsc;
458 /* N.B. dtot and dtot2 contain the total distance and the total squared
459 * distance respectively, BUT they return RMS and the scaled RMS resp.
460 */
461
462 *rmsmax=-1000;
463 *rmscmax=-1000;
464 *meanmax=-1000;
465
466 for(i=0; (i<nind-1); i++) {
467 for(j=i+1; (j<nind); j++) {
468 mean =dtot[i][j]/nframes;
469 mean2=dtot2[i][j]/nframes;
470 rms=sqrt(max(0,mean2-mean*mean));
471 rmsc=rms/mean;
472 if (mean > *meanmax) *meanmax=mean;
473 if (rms > *rmsmax ) *rmsmax =rms;
474 if (rmsc > *rmscmax) *rmscmax=rmsc;
475 meanmat[i][j]=meanmat[j][i]=mean;
476 rmsmat[i][j] =rmsmat[j][i] =rms;
477 rmscmat[i][j]=rmscmat[j][i]=rmsc;
478 }
479 }
480 }
481
482 real rms_diff(int natom,real **d,real **d_r)
483 {
484 int i,j;
485 real r,r2;
486
487 r2=0.0;
488 for(i=0; (i<natom-1); i++)
489 for(j=i+1; (j<natom); j++) {
490 r=d[i][j]-d_r[i][j];
491 r2+=r*r;
492 }
493 r2/=(natom*(natom-1))/2;
494
495 return sqrt(r2);
496 }
497
498 int gmx_rmsdist (int argc,char *argv[])
499 {
500 static char *desc[] = {
501 "g_rmsdist computes the root mean square deviation of atom distances,",
502 "which has the advantage that no fit is needed like in standard RMS",
503 "deviation as computed by g_rms.",
504 "The reference structure is taken from the structure file.",
505 "The rmsd at time t is calculated as the rms",
506 "of the differences in distance between atom-pairs in the reference",
507 "structure and the structure at time t.[PAR]",
508 "g_rmsdist can also produce matrices of the rms distances, rms distances",
509 "scaled with the mean distance and the mean distances and matrices with",
510 "NMR averaged distances (1/r^3 and 1/r^6 averaging). Finally, lists",
511 "of atom pairs with 1/r^3 and 1/r^6 averaged distance below the",
512 "maximum distance ([TT]-max[tt], which will default to 0.6 in this case)",
513 "can be generated, by default averaging over equivalent hydrogens",
514 "(all triplets of hydrogens named *[123]). Additionally a list of",
515 "equivalent atoms can be supplied ([TT]-equiv[tt]), each line containing",
516 "a set of equivalent atoms specified as residue number and name and",
517 "atom name; e.g.:[PAR]",
518 "[TT]3 SER HB1 3 SER HB2[tt][PAR]",
519 "Residue and atom names must exactly match those in the structure",
520 "file, including case. Specifying non-sequential atoms is undefined."
521
522 };
523
524 int natom,i,j,teller,gi,gj;
525 real t;
526
527 t_topology top;
528 t_atoms *atoms;
529 matrix box;
530 rvec *x;
531 FILE *fp;
532
533 int status,isize,gnr=0;
534 atom_id *index, *noe_index;
535 char *grpname;
536 real **d_r,**d,**dtot,**dtot2,**mean,**rms,**rmsc,*resnr;
537 real **dtot1_3=NULL,**dtot1_6=NULL;
538 real rmsnow,meanmax,rmsmax,rmscmax;
539 real max1_3,max1_6;
540 t_noe_gr *noe_gr=NULL;
541 t_noe **noe=NULL;
542 t_rgb rlo,rhi;
543 char buf[255];
544 bool bRMS, bScale, bMean, bNOE, bNMR3, bNMR6, bNMR;
545
546 static int nlevels=40;
547 static real scalemax=-1.0;
548 static bool bSumH=TRUE;
549 t_pargs pa[] = {
550 { "-nlevels", FALSE, etINT, {&nlevels},
551 "Discretize rms in # levels" },
552 { "-max", FALSE, etREAL, {&scalemax},
553 "Maximum level in matrices" },
554 { "-sumh", FALSE, etBOOL, {&bSumH},
555 "average distance over equivalent hydrogens" }
556 };
557 t_filenm fnm[] = {
558 { efTRX, "-f", NULL, ffREAD },
559 { efTPS, NULL, NULL, ffREAD },
560 { efNDX, NULL, NULL, ffOPTRD },
561 { efDAT, "-equiv","equiv", ffOPTRD },
562 { efXVG, NULL, "distrmsd", ffWRITE },
563 { efXPM, "-rms", "rmsdist", ffOPTWR },
564 { efXPM, "-scl", "rmsscale", ffOPTWR },
565 { efXPM, "-mean","rmsmean", ffOPTWR },
566 { efXPM, "-nmr3","nmr3", ffOPTWR },
567 { efXPM, "-nmr6","nmr6", ffOPTWR },
568 { efDAT, "-noe", "noe", ffOPTWR },
569 };
570 #define NFILE asize(fnm)
571
572 CopyRight(stderr,argv[0]);
573 parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE ,
574 NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL);
575
576 bRMS = opt2bSet("-rms", NFILE,fnm);
577 bScale= opt2bSet("-scl", NFILE,fnm);
578 bMean = opt2bSet("-mean",NFILE,fnm);
579 bNOE = opt2bSet("-noe", NFILE,fnm);
580 bNMR3 = opt2bSet("-nmr3",NFILE,fnm);
581 bNMR6 = opt2bSet("-nmr6",NFILE,fnm);
582 bNMR = bNMR3 || bNMR6 || bNOE;
583
584 /* check input */
585 if (bNOE && scalemax < 0) {
586 scalemax=0.6;
587 fprintf(stderr,"WARNING: using -noe without -max "
588 "makes no sense, setting -max to %g\n\n",scalemax);
589 }
590
591 /* get topology and index */
592 read_tps_conf(ftp2fn(efTPS,NFILE,fnm),buf,&top,&x,NULL,box,FALSE);
593 atoms=&(top.atoms);
594
595 get_index(atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&isize,&index,&grpname);
596
597 /* initialize arrays */
598 snew(d,isize);
599 snew(dtot,isize);
600 snew(dtot2,isize);
601 if (bNMR) {
602 snew(dtot1_3,isize);
603 snew(dtot1_6,isize);
604 }
605 snew(mean,isize);
606 snew(rms,isize);
607 snew(rmsc,isize);
608 snew(d_r,isize);
609 snew(resnr,isize);
610 for(i=0; (i<isize); i++) {
611 snew(d[i],isize);
612 snew(dtot[i],isize);
613 snew(dtot2[i],isize);
614 if (bNMR) {
615 snew(dtot1_3[i],isize);
616 snew(dtot1_6[i],isize);
617 }
618 snew(mean[i],isize);
619 snew(rms[i],isize);
620 snew(rmsc[i],isize);
621 snew(d_r[i],isize);
622 resnr[i]=i+1;
623 }
624
625 /*set box type*/
626 init_pbc(box);
627 calc_dist(isize,index,x,d_r);
628 sfree(x);
629
630 /*open output files*/
631 fp=xvgropen(ftp2fn(efXVG,NFILE,fnm),"RMS Deviation","Time (ps)","RMSD (nm)");
632 fprintf(fp,"@ subtitle \"of distances between %s atoms\"\n",grpname);
633
634 /*do a first step*/
635 natom=read_first_x(&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
636
637 do {
638 calc_dist_tot(isize,index,x,d,dtot,dtot2,bNMR,dtot1_3,dtot1_6);
639
640 rmsnow=rms_diff(isize,d,d_r);
641 fprintf(fp,"%g %g\n",t,rmsnow);
642 } while (read_next_x(status,&t,natom,x,box));
643 fprintf(stderr, "\n");
644
645 fclose(fp);
646 close_trj(status);
647
648 teller = nframes_read();
649 calc_rms(isize,teller,dtot,dtot2,mean,&meanmax,rms,&rmsmax,rmsc,&rmscmax);
650 fprintf(stderr,"rmsmax = %g, rmscmax = %g\n",rmsmax,rmscmax);
651
652 if (bNMR) {
653 max1_3=0;
654 max1_6=0;
655 calc_nmr(isize,teller,dtot1_3,dtot1_6,&max1_3,&max1_6);
656 }
657
658 if (scalemax > -1.0) {
659 rmsmax=scalemax;
660 rmscmax=scalemax;
661 meanmax=scalemax;
662 max1_3=scalemax;
663 max1_6=scalemax;
664 }
665
666 if (bNOE) {
667 /* make list of noe atom groups */
668 snew(noe_index, isize+1);
669 snew(noe_gr, isize);
670 gnr=analyze_noe_equivalent(opt2fn_null("-equiv",NFILE,fnm),
671 atoms, isize, index, bSumH, noe_index, noe_gr);
672 fprintf(stdout, "Found %d non-equivalent atom-groups in %d atoms\n",
673 gnr, isize);
674 /* make half matrix of of noe-group distances from atom distances */
675 snew(noe,gnr);
676 for(i=0; i<gnr; i++)
677 snew(noe[i],gnr);
678 calc_noe(isize, noe_index, dtot1_3, dtot1_6, gnr, noe);
679 }
680
681 rlo.r=1.0, rlo.g=1.0, rlo.b=1.0;
682 rhi.r=0.0, rhi.g=0.0, rhi.b=0.0;
683
684 if ( bRMS )
685 write_xpm(opt2FILE("-rms",NFILE,fnm,"w"),
686 "RMS of distance","RMS (nm)","Atom Index","Atom Index",
687 isize,isize,resnr,resnr,rms,0.0,rmsmax,rlo,rhi,&nlevels);
688
689 if ( bScale )
690 write_xpm(opt2FILE("-scl",NFILE,fnm,"w"),
691 "Relative RMS","RMS","Atom Index","Atom Index",
692 isize,isize,resnr,resnr,rmsc,0.0,rmscmax,rlo,rhi,&nlevels);
693
694 if ( bMean )
695 write_xpm(opt2FILE("-mean",NFILE,fnm,"w"),
696 "Mean Distance","Distance (nm)","Atom Index","Atom Index",
697 isize,isize,resnr,resnr,mean,0.0,meanmax,rlo,rhi,&nlevels);
698
699 if (bNMR3)
700 write_xpm(opt2FILE("-nmr3",NFILE,fnm,"w"),"1/r^3 averaged distances",
701 "Distance (nm)","Atom Index","Atom Index",
702 isize,isize,resnr,resnr,dtot1_3,0.0,max1_3,rlo,rhi,&nlevels);
703 if (bNMR6)
704 write_xpm(opt2FILE("-nmr6",NFILE,fnm,"w"),"1/r^6 averaged distances",
705 "Distance (nm)","Atom Index","Atom Index",
706 isize,isize,resnr,resnr,dtot1_6,0.0,max1_6,rlo,rhi,&nlevels);
707
708 if (bNOE)
709 write_noe(opt2FILE("-noe",NFILE,fnm,"w"), gnr, noe, noe_gr, scalemax);
710
711 do_view(ftp2fn(efXVG,NFILE,fnm),NULL);
712
713 thanx(stderr);
714 return 0;
715 }
The ultimate molecular dynamics simulation package