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