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