3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
60 static void periodic_dist(matrix box
,rvec x
[],int n
,atom_id index
[],
61 real
*rmin
,real
*rmax
,int *min_ind
)
65 real sqr_box
,r2min
,r2max
,r2
;
66 rvec shift
[NSHIFT
],d0
,d
;
68 sqr_box
= sqr(min(box
[XX
][XX
],min(box
[YY
][YY
],box
[ZZ
][ZZ
])));
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) {
76 shift
[s
][i
] = sx
*box
[XX
][i
]+sy
*box
[YY
][i
]+sz
*box
[ZZ
][i
];
84 for(j
=i
+1; j
<n
; j
++) {
85 rvec_sub(x
[index
[i
]],x
[index
[j
]],d0
);
89 for(s
=0; s
<NSHIFT
; s
++) {
90 rvec_add(d0
,shift
[s
],d
);
104 static void periodic_mindist_plot(const char *trxfn
,const char *outfn
,
105 t_topology
*top
,int ePBC
,
106 int n
,atom_id index
[],bool bSplit
,
107 const output_env_t oenv
)
110 char *leg
[5] = { "min per.","max int.","box1","box2","box3" };
115 int natoms
,ind_min
[2]={0,0},ind_mini
=0,ind_minj
=0;
116 real r
,rmin
,rmax
,rmint
,tmint
;
119 natoms
=read_first_x(oenv
,&status
,trxfn
,&t
,&x
,box
);
121 check_index(NULL
,n
,index
,NULL
,natoms
);
123 out
= xvgropen(outfn
,"Minimum distance to periodic image",
124 output_env_get_time_label(oenv
),"Distance (nm)",oenv
);
125 if (output_env_get_print_xvgr_codes(oenv
))
126 fprintf(out
,"@ subtitle \"and maximum internal distance\"\n");
127 xvgr_legend(out
,5,leg
,oenv
);
135 rm_pbc(&(top
->idef
),ePBC
,natoms
,box
,x
,x
);
137 periodic_dist(box
,x
,n
,index
,&rmin
,&rmax
,ind_min
);
141 ind_mini
= ind_min
[0];
142 ind_minj
= ind_min
[1];
144 if ( bSplit
&& !bFirst
&& abs(t
/output_env_get_time_factor(oenv
))<1e-5 )
146 fprintf(out
,"\t%g\t%6.3f %6.3f %6.3f %6.3f %6.3f\n",
147 output_env_conv_time(oenv
,t
),rmin
,rmax
,norm(box
[0]),norm(box
[1]),norm(box
[2]));
149 } while(read_next_x(oenv
,status
,&t
,natoms
,x
,box
));
154 "\nThe shortest periodic distance is %g (nm) at time %g (%s),\n"
155 "between atoms %d and %d\n",
156 rmint
,output_env_conv_time(oenv
,tmint
),output_env_get_time_unit(oenv
),
157 index
[ind_mini
]+1,index
[ind_minj
]+1);
160 static void calc_dist(real rcut
, bool bPBC
, int ePBC
, matrix box
, rvec x
[],
161 int nx1
,int nx2
, atom_id index1
[], atom_id index2
[],
163 real
*rmin
, real
*rmax
, int *nmin
, int *nmax
,
164 int *ixmin
, int *jxmin
, int *ixmax
, int *jxmax
)
170 real r2
,rmin2
,rmax2
,rcut2
;
183 /* Must init pbc every step because of pressure coupling */
185 set_pbc(&pbc
,ePBC
,box
);
198 for(j
=0; (j
< j1
); j
++) {
200 if (index2
== NULL
) {
205 for(i
=i0
; (i
< nx1
); i
++) {
209 pbc_dx(&pbc
,x
[ix
],x
[jx
],dx
);
211 rvec_sub(x
[ix
],x
[jx
],dx
);
225 } else if (r2
> rcut2
) {
246 void dist_plot(const char *fn
,const char *afile
,const char *dfile
,
247 const char *nfile
,const char *rfile
,const char *xfile
,
248 real rcut
,bool bMat
,t_atoms
*atoms
,
249 int ng
,atom_id
*index
[],int gnx
[],char *grpn
[],bool bSplit
,
250 bool bMin
, int nres
, atom_id
*residue
,bool bPBC
,int ePBC
,
251 bool bGroup
,bool bEachResEachTime
, bool bPrintResName
,
252 const output_env_t oenv
)
254 FILE *atm
,*dist
,*num
;
258 real t
,dmin
,dmax
,**mindres
=NULL
,**maxdres
=NULL
;
259 int nmin
,nmax
,status
;
261 int min1
,min2
,max1
,max2
;
267 FILE *respertime
=NULL
;
269 if ((natoms
=read_first_x(oenv
,&status
,fn
,&t
,&x0
,box
))==0)
270 gmx_fatal(FARGS
,"Could not read coordinates from statusfile\n");
272 sprintf(buf
,"%simum Distance",bMin
? "Min" : "Max");
273 dist
= xvgropen(dfile
,buf
,output_env_get_time_label(oenv
),"Distance (nm)",oenv
);
274 sprintf(buf
,"Number of Contacts %s %g nm",bMin
? "<" : ">",rcut
);
275 num
= nfile
? xvgropen(nfile
,buf
,output_env_get_time_label(oenv
),"Number",oenv
) : NULL
;
276 atm
= afile
? ffopen(afile
,"w") : NULL
;
277 trxout
= xfile
? open_trx(xfile
,"w") : NOTSET
;
282 sprintf(buf
,"Internal in %s",grpn
[0]);
284 xvgr_legend(dist
,0,leg
,oenv
);
285 if (num
) xvgr_legend(num
,0,leg
,oenv
);
288 snew(leg
,(ng
*(ng
-1))/2);
289 for(i
=j
=0; (i
<ng
-1); i
++) {
290 for(k
=i
+1; (k
<ng
); k
++,j
++) {
291 sprintf(buf
,"%s-%s",grpn
[i
],grpn
[k
]);
295 xvgr_legend(dist
,j
,leg
,oenv
);
296 if (num
) xvgr_legend(num
,j
,leg
,oenv
);
301 for(i
=0; (i
<ng
-1); i
++) {
302 sprintf(buf
,"%s-%s",grpn
[0],grpn
[i
+1]);
305 xvgr_legend(dist
,ng
-1,leg
,oenv
);
306 if (num
) xvgr_legend(num
,ng
-1,leg
,oenv
);
309 if (bEachResEachTime
)
311 sprintf(buf
,"%simum Distance",bMin
? "Min" : "Max");
312 respertime
=xvgropen(rfile
,buf
,output_env_get_time_label(oenv
),"Distance (nm)",oenv
);
313 xvgr_legend(respertime
,ng
-1,leg
,oenv
);
315 fprintf(respertime
,"# ");
316 for (j
=0; j
<nres
; j
++)
317 fprintf(respertime
,"%s%d ",*(atoms
->resinfo
[atoms
->atom
[index
[0][residue
[j
]]].resind
].name
),atoms
->atom
[index
[0][residue
[j
]]].resind
);
318 fprintf(respertime
,"\n");
326 for(i
=1; i
<ng
; i
++) {
327 snew(mindres
[i
-1], nres
);
328 snew(maxdres
[i
-1], nres
);
329 for(j
=0; j
<nres
; j
++)
331 /* maxdres[*][*] is already 0 */
336 if ( bSplit
&& !bFirst
&& abs(t
/output_env_get_time_factor(oenv
))<1e-5 ) {
337 fprintf(dist
, "&\n");
338 if (num
) fprintf(num
, "&\n");
339 if (atm
) fprintf(atm
, "&\n");
341 fprintf(dist
,"%12e",output_env_conv_time(oenv
,t
));
342 if (num
) fprintf(num
,"%12e",output_env_conv_time(oenv
,t
));
346 calc_dist(rcut
,bPBC
,ePBC
,box
,x0
,gnx
[0],gnx
[0],index
[0],index
[0],bGroup
,
347 &dmin
,&dmax
,&nmin
,&nmax
,&min1
,&min2
,&max1
,&max2
);
348 fprintf(dist
," %12e",bMin
?dmin
:dmax
);
349 if (num
) fprintf(num
," %8d",bMin
?nmin
:nmax
);
352 for(i
=0; (i
<ng
-1); i
++) {
353 for(k
=i
+1; (k
<ng
); k
++) {
354 calc_dist(rcut
,bPBC
,ePBC
,box
,x0
,gnx
[i
],gnx
[k
],index
[i
],index
[k
],
355 bGroup
,&dmin
,&dmax
,&nmin
,&nmax
,&min1
,&min2
,&max1
,&max2
);
356 fprintf(dist
," %12e",bMin
?dmin
:dmax
);
357 if (num
) fprintf(num
," %8d",bMin
?nmin
:nmax
);
363 for(i
=1; (i
<ng
); i
++) {
364 calc_dist(rcut
,bPBC
,ePBC
,box
,x0
,gnx
[0],gnx
[i
],index
[0],index
[i
],bGroup
,
365 &dmin
,&dmax
,&nmin
,&nmax
,&min1
,&min2
,&max1
,&max2
);
366 fprintf(dist
," %12e",bMin
?dmin
:dmax
);
367 if (num
) fprintf(num
," %8d",bMin
?nmin
:nmax
);
369 for(j
=0; j
<nres
; j
++) {
370 calc_dist(rcut
,bPBC
,ePBC
,box
,x0
,residue
[j
+1]-residue
[j
],gnx
[i
],
371 &(index
[0][residue
[j
]]),index
[i
],bGroup
,
372 &dmin
,&dmax
,&nmin
,&nmax
,&min1
,&min2
,&max1
,&max2
);
373 mindres
[i
-1][j
] = min(mindres
[i
-1][j
],dmin
);
374 maxdres
[i
-1][j
] = max(maxdres
[i
-1][j
],dmax
);
382 if ( bMin
?min1
:max1
!= -1 )
384 fprintf(atm
,"%12e %12d %12d\n",
385 output_env_conv_time(oenv
,t
),1+(bMin
? min1
: max1
),
386 1+(bMin
? min2
: max2
));
389 oindex
[0]=bMin
?min1
:max1
;
390 oindex
[1]=bMin
?min2
:max2
;
391 write_trx(trxout
,2,oindex
,atoms
,i
,t
,box
,x0
,NULL
,NULL
);
394 /*dmin should be minimum distance for residue and group*/
395 if (bEachResEachTime
)
397 fprintf(respertime
,"%12e",t
);
399 for(j
=0; j
<nres
; j
++)
401 fprintf(respertime
, " %7g", bMin
? mindres
[i
-1][j
] : maxdres
[i
-1][j
]);
402 /*reset distances for next time point*/
406 fprintf(respertime
, "\n");
408 } while (read_next_x(oenv
,status
,&t
,natoms
,x0
,box
));
412 if (num
) ffclose(num
);
413 if (atm
) ffclose(atm
);
414 if (trxout
>=0) close_xtc(trxout
);
416 if(nres
&& !bEachResEachTime
) {
419 sprintf(buf
,"%simum Distance",bMin
? "Min" : "Max");
420 res
=xvgropen(rfile
,buf
,"Residue (#)","Distance (nm)",oenv
);
421 xvgr_legend(res
,ng
-1,leg
,oenv
);
422 for(j
=0; j
<nres
; j
++) {
423 fprintf(res
, "%4d", j
+1);
424 for(i
=1; i
<ng
; i
++) {
425 fprintf(res
, " %7g", bMin
? mindres
[i
-1][j
] : maxdres
[i
-1][j
]);
434 int find_residues(t_atoms
*atoms
, int n
, atom_id index
[], atom_id
**resindex
)
437 int nres
=0,resnr
, presnr
;
440 /* build index of first atom numbers for each residue */
442 snew(residx
, atoms
->nres
+1);
444 resnr
= atoms
->atom
[index
[i
]].resind
;
445 if (resnr
!= presnr
) {
451 if (debug
) printf("Found %d residues out of %d (%d/%d atoms)\n",
452 nres
, atoms
->nres
, atoms
->nr
, n
);
453 srenew(residx
, nres
+1);
454 /* mark end of last residue */
460 void dump_res(FILE *out
, int nres
, atom_id
*resindex
, int n
, atom_id index
[])
464 for(i
=0; i
<nres
-1; i
++) {
465 fprintf(out
,"Res %d (%d):", i
, resindex
[i
+1]-resindex
[i
]);
466 for(j
=resindex
[i
]; j
<resindex
[i
+1]; j
++)
467 fprintf(out
," %d(%d)", j
, index
[j
]);
472 int gmx_mindist(int argc
,char *argv
[])
474 const char *desc
[] = {
475 "g_mindist computes the distance between one group and a number of",
476 "other groups. Both the minimum distance",
477 "(between any pair of atoms from the respective groups)",
478 "and the number of contacts within a given",
479 "distance are written to two separate output files.",
480 "With the [TT]-group[tt] option a contact of an atom an other group",
481 "with multiple atoms in the first group is counted as one contact",
482 "instead of as multiple contacts.",
483 "With [TT]-or[tt], minimum distances to each residue in the first",
484 "group are determined and plotted as a function of reisdue number.[PAR]",
485 "With option [TT]-pi[tt] the minimum distance of a group to its",
486 "periodic image is plotted. This is useful for checking if a protein",
487 "has seen its periodic image during a simulation. Only one shift in",
488 "each direction is considered, giving a total of 26 shifts.",
489 "It also plots the maximum distance within the group and the lengths",
490 "of the three box vectors.[PAR]",
491 "Other programs that calculate distances are [TT]g_dist[tt]",
492 "and [TT]g_bond[tt]."
494 const char *bugs
[] = {
495 "The [TT]-pi[tt] option is very slow."
498 static bool bMat
=FALSE
,bPI
=FALSE
,bSplit
=FALSE
,bMax
=FALSE
,bPBC
=TRUE
;
499 static bool bGroup
=FALSE
;
500 static real rcutoff
=0.6;
502 static bool bEachResEachTime
=FALSE
,bPrintResName
=FALSE
;
504 { "-matrix", FALSE
, etBOOL
, {&bMat
},
505 "Calculate half a matrix of group-group distances" },
506 { "-max", FALSE
, etBOOL
, {&bMax
},
507 "Calculate *maximum* distance instead of minimum" },
508 { "-d", FALSE
, etREAL
, {&rcutoff
},
509 "Distance for contacts" },
510 { "-group", FALSE
, etBOOL
, {&bGroup
},
511 "Count contacts with multiple atoms in the first group as one" },
512 { "-pi", FALSE
, etBOOL
, {&bPI
},
513 "Calculate minimum distance with periodic images" },
514 { "-split", FALSE
, etBOOL
, {&bSplit
},
515 "Split graph where time is zero" },
516 { "-ng", FALSE
, etINT
, {&ng
},
517 "Number of secondary groups to compute distance to a central group" },
518 { "-pbc", FALSE
, etBOOL
, {&bPBC
},
519 "Take periodic boundary conditions into account" },
520 { "-respertime", FALSE
, etBOOL
, {&bEachResEachTime
},
521 "When writing per-residue distances, write distance for each time point" },
522 { "-printresname", FALSE
, etBOOL
, {&bPrintResName
},
523 "Write residue names" }
526 t_topology
*top
=NULL
;
536 const char *trxfnm
,*tpsfnm
,*ndxfnm
,*distfnm
,*numfnm
,*atmfnm
,*oxfnm
,*resfnm
;
539 atom_id
**index
, *residues
=NULL
;
541 { efTRX
, "-f", NULL
, ffREAD
},
542 { efTPS
, NULL
, NULL
, ffOPTRD
},
543 { efNDX
, NULL
, NULL
, ffOPTRD
},
544 { efXVG
, "-od","mindist", ffWRITE
},
545 { efXVG
, "-on","numcont", ffOPTWR
},
546 { efOUT
, "-o", "atm-pair", ffOPTWR
},
547 { efTRO
, "-ox","mindist", ffOPTWR
},
548 { efXVG
, "-or","mindistres", ffOPTWR
}
550 #define NFILE asize(fnm)
552 CopyRight(stderr
,argv
[0]);
553 parse_common_args(&argc
,argv
,
554 PCA_CAN_VIEW
| PCA_CAN_TIME
| PCA_TIME_UNIT
| PCA_BE_NICE
,
555 NFILE
,fnm
,asize(pa
),pa
,asize(desc
),desc
,0,NULL
, &oenv
);
557 trxfnm
= ftp2fn(efTRX
,NFILE
,fnm
);
558 ndxfnm
= ftp2fn_null(efNDX
,NFILE
,fnm
);
559 distfnm
= opt2fn("-od",NFILE
,fnm
);
560 numfnm
= opt2fn_null("-on",NFILE
,fnm
);
561 atmfnm
= ftp2fn_null(efOUT
,NFILE
,fnm
);
562 oxfnm
= opt2fn_null("-ox",NFILE
,fnm
);
563 resfnm
= opt2fn_null("-or",NFILE
,fnm
);
564 if (bPI
|| resfnm
!= NULL
) {
565 /* We need a tps file */
566 tpsfnm
= ftp2fn(efTPS
,NFILE
,fnm
);
568 tpsfnm
= ftp2fn_null(efTPS
,NFILE
,fnm
);
571 if (!tpsfnm
&& !ndxfnm
)
572 gmx_fatal(FARGS
,"You have to specify either the index file or a tpr file");
576 fprintf(stderr
,"Choose a group for distance calculation\n");
585 if (tpsfnm
|| resfnm
|| !ndxfnm
) {
587 bTop
= read_tps_conf(tpsfnm
,title
,top
,&ePBC
,&x
,NULL
,box
,FALSE
);
589 printf("\nWARNING: Without a run input file a trajectory with broken molecules will not give the correct periodic image distance\n\n");
591 get_index(top
? &(top
->atoms
) : NULL
,ndxfnm
,ng
,gnx
,index
,grpname
);
593 if (bMat
&& (ng
== 1)) {
595 printf("Special case: making distance matrix between all atoms in group %s\n",
600 for(i
=1; (i
<ng
); i
++) {
602 grpname
[i
] = grpname
[0];
604 index
[i
][0] = index
[0][i
];
610 nres
=find_residues(top
? &(top
->atoms
) : NULL
,
611 gnx
[0], index
[0], &residues
);
612 if (debug
) dump_res(debug
, nres
, residues
, gnx
[0], index
[0]);
616 periodic_mindist_plot(trxfnm
,distfnm
,top
,ePBC
,gnx
[0],index
[0],bSplit
,oenv
);
618 dist_plot(trxfnm
,atmfnm
,distfnm
,numfnm
,resfnm
,oxfnm
,
619 rcutoff
,bMat
,top
? &(top
->atoms
) : NULL
,
620 ng
,index
,gnx
,grpname
,bSplit
,!bMax
, nres
, residues
,bPBC
,ePBC
,
621 bGroup
,bEachResEachTime
,bPrintResName
,oenv
);
624 do_view(oenv
,distfnm
,"-nxy");
626 do_view(oenv
,numfnm
,"-nxy");