?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Fixes #1183 PBC bug in g_mindist
[gromacs/AngularHB.git] / src / tools / gmx_mindist.c
blobdc9ae1f74f52055bd7a76c0154029148c9d0e111
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 <stdlib.h>
41
42 #include "sysstuff.h"
43 #include "string.h"
44 #include "typedefs.h"
45 #include "smalloc.h"
46 #include "macros.h"
47 #include "vec.h"
48 #include "xvgr.h"
49 #include "pbc.h"
50 #include "copyrite.h"
51 #include "futil.h"
52 #include "statutil.h"
53 #include "index.h"
54 #include "tpxio.h"
55 #include "rmpbc.h"
56 #include "xtcio.h"
57 #include "gmx_ana.h"
58
59
60 static void periodic_dist(matrix box,rvec x[],int n,atom_id index[],
61 real *rmin,real *rmax,int *min_ind)
62 {
63 #define NSHIFT 26
64 int sx,sy,sz,i,j,s;
65 real sqr_box,r2min,r2max,r2;
66 rvec shift[NSHIFT],d0,d;
67
68 sqr_box = sqr(min(norm(box[XX]),min(norm(box[YY]),norm(box[ZZ]))));
69
70 s = 0;
71 for(sz=-1; sz<=1; sz++)
72 for(sy=-1; sy<=1; sy++)
73 for(sx=-1; sx<=1; sx++)
74 if (sx!=0 || sy!=0 || sz!=0) {
75 for(i=0; i<DIM; i++)
76 shift[s][i] = sx*box[XX][i]+sy*box[YY][i]+sz*box[ZZ][i];
77 s++;
78 }
79
80 r2min = sqr_box;
81 r2max = 0;
82
83 for(i=0; i<n; i++)
84 for(j=i+1; j<n; j++) {
85 rvec_sub(x[index[i]],x[index[j]],d0);
86 r2 = norm2(d0);
87 if (r2 > r2max)
88 r2max = r2;
89 for(s=0; s<NSHIFT; s++) {
90 rvec_add(d0,shift[s],d);
91 r2 = norm2(d);
92 if (r2 < r2min) {
93 r2min = r2;
94 min_ind[0] = i;
95 min_ind[1] = j;
96 }
97 }
98 }
99
100 *rmin = sqrt(r2min);
101 *rmax = sqrt(r2max);
102 }
103
104 static void periodic_mindist_plot(const char *trxfn,const char *outfn,
105 t_topology *top,int ePBC,
106 int n,atom_id index[],gmx_bool bSplit,
107 const output_env_t oenv)
108 {
109 FILE *out;
110 const char *leg[5] = { "min per.","max int.","box1","box2","box3" };
111 t_trxstatus *status;
112 real t;
113 rvec *x;
114 matrix box;
115 int natoms,ind_min[2]={0,0},ind_mini=0,ind_minj=0;
116 real r,rmin,rmax,rmint,tmint;
117 gmx_bool bFirst;
118 gmx_rmpbc_t gpbc=NULL;
119
120 natoms=read_first_x(oenv,&status,trxfn,&t,&x,box);
121
122 check_index(NULL,n,index,NULL,natoms);
123
124 out = xvgropen(outfn,"Minimum distance to periodic image",
125 output_env_get_time_label(oenv),"Distance (nm)",oenv);
126 if (output_env_get_print_xvgr_codes(oenv))
127 fprintf(out,"@ subtitle \"and maximum internal distance\"\n");
128 xvgr_legend(out,5,leg,oenv);
129
130 rmint = box[XX][XX];
131 tmint = 0;
132
133 if (NULL != top)
134 gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box);
135
136 bFirst=TRUE;
137 do {
138 if (NULL != top)
139 gmx_rmpbc(gpbc,natoms,box,x);
140
141 periodic_dist(box,x,n,index,&rmin,&rmax,ind_min);
142 if (rmin < rmint) {
143 rmint = rmin;
144 tmint = t;
145 ind_mini = ind_min[0];
146 ind_minj = ind_min[1];
147 }
148 if ( bSplit && !bFirst && abs(t/output_env_get_time_factor(oenv))<1e-5 )
149 fprintf(out, "&\n");
150 fprintf(out,"\t%g\t%6.3f %6.3f %6.3f %6.3f %6.3f\n",
151 output_env_conv_time(oenv,t),rmin,rmax,norm(box[0]),norm(box[1]),norm(box[2]));
152 bFirst=FALSE;
153 } while(read_next_x(oenv,status,&t,natoms,x,box));
154
155 if (NULL != top)
156 gmx_rmpbc_done(gpbc);
157
158 ffclose(out);
159
160 fprintf(stdout,
161 "\nThe shortest periodic distance is %g (nm) at time %g (%s),\n"
162 "between atoms %d and %d\n",
163 rmint,output_env_conv_time(oenv,tmint),output_env_get_time_unit(oenv),
164 index[ind_mini]+1,index[ind_minj]+1);
165 }
166
167 static void calc_dist(real rcut, gmx_bool bPBC, int ePBC, matrix box, rvec x[],
168 int nx1,int nx2, atom_id index1[], atom_id index2[],
169 gmx_bool bGroup,
170 real *rmin, real *rmax, int *nmin, int *nmax,
171 int *ixmin, int *jxmin, int *ixmax, int *jxmax)
172 {
173 int i,j,i0=0,j1;
174 int ix,jx;
175 atom_id *index3;
176 rvec dx;
177 real r2,rmin2,rmax2,rcut2;
178 t_pbc pbc;
179 int nmin_j,nmax_j;
180
181 *ixmin = -1;
182 *jxmin = -1;
183 *ixmax = -1;
184 *jxmax = -1;
185 *nmin = 0;
186 *nmax = 0;
187
188 rcut2=sqr(rcut);
189
190 /* Must init pbc every step because of pressure coupling */
191 if (bPBC)
192 set_pbc(&pbc,ePBC,box);
193 if (index2) {
194 i0=0;
195 j1=nx2;
196 index3=index2;
197 } else {
198 j1=nx1;
199 index3=index1;
200 }
201
202 rmin2=1e12;
203 rmax2=-1e12;
204
205 for(j=0; (j < j1); j++) {
206 jx = index3[j];
207 if (index2 == NULL) {
208 i0 = j + 1;
209 }
210 nmin_j = 0;
211 nmax_j = 0;
212 for(i=i0; (i < nx1); i++) {
213 ix = index1[i];
214 if (ix != jx) {
215 if (bPBC)
216 pbc_dx(&pbc,x[ix],x[jx],dx);
217 else
218 rvec_sub(x[ix],x[jx],dx);
219 r2=iprod(dx,dx);
220 if (r2 < rmin2) {
221 rmin2=r2;
222 *ixmin=ix;
223 *jxmin=jx;
224 }
225 if (r2 > rmax2) {
226 rmax2=r2;
227 *ixmax=ix;
228 *jxmax=jx;
229 }
230 if (r2 <= rcut2) {
231 nmin_j++;
232 } else if (r2 > rcut2) {
233 nmax_j++;
234 }
235 }
236 }
237 if (bGroup) {
238 if (nmin_j > 0) {
239 (*nmin)++;
240 }
241 if (nmax_j > 0) {
242 (*nmax)++;
243 }
244 } else {
245 *nmin += nmin_j;
246 *nmax += nmax_j;
247 }
248 }
249 *rmin = sqrt(rmin2);
250 *rmax = sqrt(rmax2);
251 }
252
253 void dist_plot(const char *fn,const char *afile,const char *dfile,
254 const char *nfile,const char *rfile,const char *xfile,
255 real rcut,gmx_bool bMat,t_atoms *atoms,
256 int ng,atom_id *index[],int gnx[],char *grpn[],gmx_bool bSplit,
257 gmx_bool bMin, int nres, atom_id *residue,gmx_bool bPBC,int ePBC,
258 gmx_bool bGroup,gmx_bool bEachResEachTime, gmx_bool bPrintResName,
259 const output_env_t oenv)
260 {
261 FILE *atm,*dist,*num;
262 t_trxstatus *trxout;
263 char buf[256];
264 char **leg;
265 real t,dmin,dmax,**mindres=NULL,**maxdres=NULL;
266 int nmin,nmax;
267 t_trxstatus *status;
268 int i=-1,j,k,natoms;
269 int min1,min2,max1,max2,min1r,min2r,max1r,max2r;
270 atom_id oindex[2];
271 rvec *x0;
272 matrix box;
273 t_trxframe frout;
274 gmx_bool bFirst;
275 FILE *respertime=NULL;
276
277 if ((natoms=read_first_x(oenv,&status,fn,&t,&x0,box))==0)
278 gmx_fatal(FARGS,"Could not read coordinates from statusfile\n");
279
280 sprintf(buf,"%simum Distance",bMin ? "Min" : "Max");
281 dist= xvgropen(dfile,buf,output_env_get_time_label(oenv),"Distance (nm)",oenv);
282 sprintf(buf,"Number of Contacts %s %g nm",bMin ? "<" : ">",rcut);
283 num = nfile ? xvgropen(nfile,buf,output_env_get_time_label(oenv),"Number",oenv) : NULL;
284 atm = afile ? ffopen(afile,"w") : NULL;
285 trxout = xfile ? open_trx(xfile,"w") : NULL;
286
287 if (bMat) {
288 if (ng == 1) {
289 snew(leg,1);
290 sprintf(buf,"Internal in %s",grpn[0]);
291 leg[0]=strdup(buf);
292 xvgr_legend(dist,0,(const char**)leg,oenv);
293 if (num) xvgr_legend(num,0,(const char**)leg,oenv);
294 }
295 else {
296 snew(leg,(ng*(ng-1))/2);
297 for(i=j=0; (i<ng-1); i++) {
298 for(k=i+1; (k<ng); k++,j++) {
299 sprintf(buf,"%s-%s",grpn[i],grpn[k]);
300 leg[j]=strdup(buf);
301 }
302 }
303 xvgr_legend(dist,j,(const char**)leg,oenv);
304 if (num) xvgr_legend(num,j,(const char**)leg,oenv);
305 }
306 }
307 else {
308 snew(leg,ng-1);
309 for(i=0; (i<ng-1); i++) {
310 sprintf(buf,"%s-%s",grpn[0],grpn[i+1]);
311 leg[i]=strdup(buf);
312 }
313 xvgr_legend(dist,ng-1,(const char**)leg,oenv);
314 if (num) xvgr_legend(num,ng-1,(const char**)leg,oenv);
315 }
316
317 if (bEachResEachTime)
318 {
319 sprintf(buf,"%simum Distance",bMin ? "Min" : "Max");
320 respertime=xvgropen(rfile,buf,output_env_get_time_label(oenv),"Distance (nm)",oenv);
321 xvgr_legend(respertime,ng-1,(const char**)leg,oenv);
322 if (bPrintResName)
323 fprintf(respertime,"# ");
324 for (j=0; j<nres; j++)
325 fprintf(respertime,"%s%d ",*(atoms->resinfo[atoms->atom[index[0][residue[j]]].resind].name),atoms->atom[index[0][residue[j]]].resind);
326 fprintf(respertime,"\n");
327
328 }
329
330 j=0;
331 if (nres) {
332 snew(mindres, ng-1);
333 snew(maxdres, ng-1);
334 for(i=1; i<ng; i++) {
335 snew(mindres[i-1], nres);
336 snew(maxdres[i-1], nres);
337 for(j=0; j<nres; j++)
338 mindres[i-1][j]=1e6;
339 /* maxdres[*][*] is already 0 */
340 }
341 }
342 bFirst=TRUE;
343 do {
344 if ( bSplit && !bFirst && abs(t/output_env_get_time_factor(oenv))<1e-5 ) {
345 fprintf(dist, "&\n");
346 if (num) fprintf(num, "&\n");
347 if (atm) fprintf(atm, "&\n");
348 }
349 fprintf(dist,"%12e",output_env_conv_time(oenv,t));
350 if (num) fprintf(num,"%12e",output_env_conv_time(oenv,t));
351
352 if (bMat) {
353 if (ng == 1) {
354 calc_dist(rcut,bPBC,ePBC,box,x0,gnx[0],gnx[0],index[0],index[0],bGroup,
355 &dmin,&dmax,&nmin,&nmax,&min1,&min2,&max1,&max2);
356 fprintf(dist," %12e",bMin?dmin:dmax);
357 if (num) fprintf(num," %8d",bMin?nmin:nmax);
358 }
359 else {
360 for(i=0; (i<ng-1); i++) {
361 for(k=i+1; (k<ng); k++) {
362 calc_dist(rcut,bPBC,ePBC,box,x0,gnx[i],gnx[k],index[i],index[k],
363 bGroup,&dmin,&dmax,&nmin,&nmax,&min1,&min2,&max1,&max2);
364 fprintf(dist," %12e",bMin?dmin:dmax);
365 if (num) fprintf(num," %8d",bMin?nmin:nmax);
366 }
367 }
368 }
369 }
370 else {
371 for(i=1; (i<ng); i++) {
372 calc_dist(rcut,bPBC,ePBC,box,x0,gnx[0],gnx[i],index[0],index[i],bGroup,
373 &dmin,&dmax,&nmin,&nmax,&min1,&min2,&max1,&max2);
374 fprintf(dist," %12e",bMin?dmin:dmax);
375 if (num) fprintf(num," %8d",bMin?nmin:nmax);
376 if (nres) {
377 for(j=0; j<nres; j++) {
378 calc_dist(rcut,bPBC,ePBC,box,x0,residue[j+1]-residue[j],gnx[i],
379 &(index[0][residue[j]]),index[i],bGroup,
380 &dmin,&dmax,&nmin,&nmax,&min1r,&min2r,&max1r,&max2r);
381 mindres[i-1][j] = min(mindres[i-1][j],dmin);
382 maxdres[i-1][j] = max(maxdres[i-1][j],dmax);
383 }
384 }
385 }
386 }
387 fprintf(dist,"\n");
388 if (num)
389 fprintf(num,"\n");
390 if ( (bMin?min1:max1) != -1 )
391 if (atm)
392 fprintf(atm,"%12e %12d %12d\n",
393 output_env_conv_time(oenv,t),1+(bMin ? min1 : max1),
394 1+(bMin ? min2 : max2));
395
396 if (trxout) {
397 oindex[0]=bMin?min1:max1;
398 oindex[1]=bMin?min2:max2;
399 write_trx(trxout,2,oindex,atoms,i,t,box,x0,NULL,NULL);
400 }
401 bFirst=FALSE;
402 /*dmin should be minimum distance for residue and group*/
403 if (bEachResEachTime)
404 {
405 fprintf(respertime,"%12e",t);
406 for(i=1; i<ng; i++)
407 for(j=0; j<nres; j++)
408 {
409 fprintf(respertime, " %7g", bMin ? mindres[i-1][j] : maxdres[i-1][j]);
410 /*reset distances for next time point*/
411 mindres[i-1][j]=1e6;
412 maxdres[i-1][j]=0;
413 }
414 fprintf(respertime, "\n");
415 }
416 } while (read_next_x(oenv,status,&t,natoms,x0,box));
417
418 close_trj(status);
419 ffclose(dist);
420 if (num) ffclose(num);
421 if (atm) ffclose(atm);
422 if (trxout) close_trx(trxout);
423
424 if(nres && !bEachResEachTime) {
425 FILE *res;
426
427 sprintf(buf,"%simum Distance",bMin ? "Min" : "Max");
428 res=xvgropen(rfile,buf,"Residue (#)","Distance (nm)",oenv);
429 xvgr_legend(res,ng-1,(const char**)leg,oenv);
430 for(j=0; j<nres; j++) {
431 fprintf(res, "%4d", j+1);
432 for(i=1; i<ng; i++) {
433 fprintf(res, " %7g", bMin ? mindres[i-1][j] : maxdres[i-1][j]);
434 }
435 fprintf(res, "\n");
436 }
437 }
438
439 sfree(x0);
440 }
441
442 int find_residues(t_atoms *atoms, int n, atom_id index[], atom_id **resindex)
443 {
444 int i;
445 int nres=0,resnr, presnr;
446 int *residx;
447
448 /* build index of first atom numbers for each residue */
449 presnr = NOTSET;
450 snew(residx, atoms->nres+1);
451 for(i=0; i<n; i++) {
452 resnr = atoms->atom[index[i]].resind;
453 if (resnr != presnr) {
454 residx[nres]=i;
455 nres++;
456 presnr=resnr;
457 }
458 }
459 if (debug) printf("Found %d residues out of %d (%d/%d atoms)\n",
460 nres, atoms->nres, atoms->nr, n);
461 srenew(residx, nres+1);
462 /* mark end of last residue */
463 residx[nres] = n;
464 *resindex = residx;
465 return nres;
466 }
467
468 void dump_res(FILE *out, int nres, atom_id *resindex, int n, atom_id index[])
469 {
470 int i,j;
471
472 for(i=0; i<nres-1; i++) {
473 fprintf(out,"Res %d (%d):", i, resindex[i+1]-resindex[i]);
474 for(j=resindex[i]; j<resindex[i+1]; j++)
475 fprintf(out," %d(%d)", j, index[j]);
476 fprintf(out,"\n");
477 }
478 }
479
480 int gmx_mindist(int argc,char *argv[])
481 {
482 const char *desc[] = {
483 "[TT]g_mindist[tt] computes the distance between one group and a number of",
484 "other groups. Both the minimum distance",
485 "(between any pair of atoms from the respective groups)",
486 "and the number of contacts within a given",
487 "distance are written to two separate output files.",
488 "With the [TT]-group[tt] option a contact of an atom in another group",
489 "with multiple atoms in the first group is counted as one contact",
490 "instead of as multiple contacts.",
491 "With [TT]-or[tt], minimum distances to each residue in the first",
492 "group are determined and plotted as a function of residue number.[PAR]",
493 "With option [TT]-pi[tt] the minimum distance of a group to its",
494 "periodic image is plotted. This is useful for checking if a protein",
495 "has seen its periodic image during a simulation. Only one shift in",
496 "each direction is considered, giving a total of 26 shifts.",
497 "It also plots the maximum distance within the group and the lengths",
498 "of the three box vectors.[PAR]",
499 "Other programs that calculate distances are [TT]g_dist[tt]",
500 "and [TT]g_bond[tt]."
501 };
502 const char *bugs[] = {
503 "The [TT]-pi[tt] option is very slow."
504 };
505
506 static gmx_bool bMat=FALSE,bPI=FALSE,bSplit=FALSE,bMax=FALSE,bPBC=TRUE;
507 static gmx_bool bGroup=FALSE;
508 static real rcutoff=0.6;
509 static int ng=1;
510 static gmx_bool bEachResEachTime=FALSE,bPrintResName=FALSE;
511 t_pargs pa[] = {
512 { "-matrix", FALSE, etBOOL, {&bMat},
513 "Calculate half a matrix of group-group distances" },
514 { "-max", FALSE, etBOOL, {&bMax},
515 "Calculate *maximum* distance instead of minimum" },
516 { "-d", FALSE, etREAL, {&rcutoff},
517 "Distance for contacts" },
518 { "-group", FALSE, etBOOL, {&bGroup},
519 "Count contacts with multiple atoms in the first group as one" },
520 { "-pi", FALSE, etBOOL, {&bPI},
521 "Calculate minimum distance with periodic images" },
522 { "-split", FALSE, etBOOL, {&bSplit},
523 "Split graph where time is zero" },
524 { "-ng", FALSE, etINT, {&ng},
525 "Number of secondary groups to compute distance to a central group" },
526 { "-pbc", FALSE, etBOOL, {&bPBC},
527 "Take periodic boundary conditions into account" },
528 { "-respertime", FALSE, etBOOL, {&bEachResEachTime},
529 "When writing per-residue distances, write distance for each time point" },
530 { "-printresname", FALSE, etBOOL, {&bPrintResName},
531 "Write residue names" }
532 };
533 output_env_t oenv;
534 t_topology *top=NULL;
535 int ePBC=-1;
536 char title[256];
537 real t;
538 rvec *x;
539 matrix box;
540 gmx_bool bTop=FALSE;
541
542 FILE *atm;
543 int i,j,nres=0;
544 const char *trxfnm,*tpsfnm,*ndxfnm,*distfnm,*numfnm,*atmfnm,*oxfnm,*resfnm;
545 char **grpname;
546 int *gnx;
547 atom_id **index, *residues=NULL;
548 t_filenm fnm[] = {
549 { efTRX, "-f", NULL, ffREAD },
550 { efTPS, NULL, NULL, ffOPTRD },
551 { efNDX, NULL, NULL, ffOPTRD },
552 { efXVG, "-od","mindist", ffWRITE },
553 { efXVG, "-on","numcont", ffOPTWR },
554 { efOUT, "-o", "atm-pair", ffOPTWR },
555 { efTRO, "-ox","mindist", ffOPTWR },
556 { efXVG, "-or","mindistres", ffOPTWR }
557 };
558 #define NFILE asize(fnm)
559
560 CopyRight(stderr,argv[0]);
561 parse_common_args(&argc,argv,
562 PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,
563 NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL, &oenv);
564
565 trxfnm = ftp2fn(efTRX,NFILE,fnm);
566 ndxfnm = ftp2fn_null(efNDX,NFILE,fnm);
567 distfnm= opt2fn("-od",NFILE,fnm);
568 numfnm = opt2fn_null("-on",NFILE,fnm);
569 atmfnm = ftp2fn_null(efOUT,NFILE,fnm);
570 oxfnm = opt2fn_null("-ox",NFILE,fnm);
571 resfnm = opt2fn_null("-or",NFILE,fnm);
572 if (bPI || resfnm != NULL) {
573 /* We need a tps file */
574 tpsfnm = ftp2fn(efTPS,NFILE,fnm);
575 } else {
576 tpsfnm = ftp2fn_null(efTPS,NFILE,fnm);
577 }
578
579 if (!tpsfnm && !ndxfnm)
580 gmx_fatal(FARGS,"You have to specify either the index file or a tpr file");
581
582 if (bPI) {
583 ng = 1;
584 fprintf(stderr,"Choose a group for distance calculation\n");
585 }
586 else if (!bMat)
587 ng++;
588
589 snew(gnx,ng);
590 snew(index,ng);
591 snew(grpname,ng);
592
593 if (tpsfnm || resfnm || !ndxfnm) {
594 snew(top,1);
595 bTop = read_tps_conf(tpsfnm,title,top,&ePBC,&x,NULL,box,FALSE);
596 if (bPI && !bTop)
597 printf("\nWARNING: Without a run input file a trajectory with broken molecules will not give the correct periodic image distance\n\n");
598 }
599 get_index(top ? &(top->atoms) : NULL,ndxfnm,ng,gnx,index,grpname);
600
601 if (bMat && (ng == 1)) {
602 ng = gnx[0];
603 printf("Special case: making distance matrix between all atoms in group %s\n",
604 grpname[0]);
605 srenew(gnx,ng);
606 srenew(index,ng);
607 srenew(grpname,ng);
608 for(i=1; (i<ng); i++) {
609 gnx[i] = 1;
610 grpname[i] = grpname[0];
611 snew(index[i],1);
612 index[i][0] = index[0][i];
613 }
614 gnx[0] = 1;
615 }
616
617 if (resfnm) {
618 nres=find_residues(top ? &(top->atoms) : NULL,
619 gnx[0], index[0], &residues);
620 if (debug) dump_res(debug, nres, residues, gnx[0], index[0]);
621 }
622
623 if (bPI) {
624 periodic_mindist_plot(trxfnm,distfnm,top,ePBC,gnx[0],index[0],bSplit,oenv);
625 } else {
626 dist_plot(trxfnm,atmfnm,distfnm,numfnm,resfnm,oxfnm,
627 rcutoff,bMat,top ? &(top->atoms) : NULL,
628 ng,index,gnx,grpname,bSplit,!bMax, nres, residues,bPBC,ePBC,
629 bGroup,bEachResEachTime,bPrintResName,oenv);
630 }
631
632 do_view(oenv,distfnm,"-nxy");
633 if (!bPI)
634 do_view(oenv,numfnm,"-nxy");
635
636 thanx(stderr);
637
638 return 0;
639 }
640