Fixed make_edi.c
[gromacs/rigid-bodies.git] / src / tools / gmx_mindist.c
blob98d4966598b6d6bbe5007fd068a34efd1ab16cd0
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.
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
32 * And Hey:
33 * Green Red Orange Magenta Azure Cyan Skyblue
35 #ifdef HAVE_CONFIG_H
36 #include <config.h>
37 #endif
39 #include <math.h>
40 #include <stdlib.h>
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"
60 static void periodic_dist(matrix box,rvec x[],int n,atom_id index[],
61 real *rmin,real *rmax,int *min_ind)
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;
68 sqr_box = sqr(min(box[XX][XX],min(box[YY][YY],box[ZZ][ZZ])));
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++;
80 r2min = sqr_box;
81 r2max = 0;
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;
100 *rmin = sqrt(r2min);
101 *rmax = sqrt(r2max);
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)
109 FILE *out;
110 char *leg[5] = { "min per.","max int.","box1","box2","box3" };
111 int 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 bool bFirst;
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);
129 rmint = box[XX][XX];
130 tmint = 0;
132 bFirst=TRUE;
133 do {
134 if (top) {
135 rm_pbc(&(top->idef),ePBC,natoms,box,x,x);
137 periodic_dist(box,x,n,index,&rmin,&rmax,ind_min);
138 if (rmin < rmint) {
139 rmint = rmin;
140 tmint = t;
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 )
145 fprintf(out, "&\n");
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]));
148 bFirst=FALSE;
149 } while(read_next_x(oenv,status,&t,natoms,x,box));
151 ffclose(out);
153 fprintf(stdout,
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[],
162 bool bGroup,
163 real *rmin, real *rmax, int *nmin, int *nmax,
164 int *ixmin, int *jxmin, int *ixmax, int *jxmax)
166 int i,j,i0=0,j1;
167 int ix,jx;
168 atom_id *index3;
169 rvec dx;
170 real r2,rmin2,rmax2,rcut2;
171 t_pbc pbc;
172 int nmin_j,nmax_j;
174 *ixmin = -1;
175 *jxmin = -1;
176 *ixmax = -1;
177 *jxmax = -1;
178 *nmin = 0;
179 *nmax = 0;
181 rcut2=sqr(rcut);
183 /* Must init pbc every step because of pressure coupling */
184 if (bPBC)
185 set_pbc(&pbc,ePBC,box);
186 if (index2) {
187 i0=0;
188 j1=nx2;
189 index3=index2;
190 } else {
191 j1=nx1;
192 index3=index1;
195 rmin2=1e12;
196 rmax2=-1e12;
198 for(j=0; (j < j1); j++) {
199 jx = index3[j];
200 if (index2 == NULL) {
201 i0 = j + 1;
203 nmin_j = 0;
204 nmax_j = 0;
205 for(i=i0; (i < nx1); i++) {
206 ix = index1[i];
207 if (ix != jx) {
208 if (bPBC)
209 pbc_dx(&pbc,x[ix],x[jx],dx);
210 else
211 rvec_sub(x[ix],x[jx],dx);
212 r2=iprod(dx,dx);
213 if (r2 < rmin2) {
214 rmin2=r2;
215 *ixmin=ix;
216 *jxmin=jx;
218 if (r2 > rmax2) {
219 rmax2=r2;
220 *ixmax=ix;
221 *jxmax=jx;
223 if (r2 <= rcut2) {
224 nmin_j++;
225 } else if (r2 > rcut2) {
226 nmax_j++;
230 if (bGroup) {
231 if (nmin_j > 0) {
232 (*nmin)++;
234 if (nmax_j > 0) {
235 (*nmax)++;
237 } else {
238 *nmin += nmin_j;
239 *nmax += nmax_j;
242 *rmin = sqrt(rmin2);
243 *rmax = sqrt(rmax2);
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;
255 int trxout;
256 char buf[256];
257 char **leg;
258 real t,dmin,dmax,**mindres=NULL,**maxdres=NULL;
259 int nmin,nmax,status;
260 int i=-1,j,k,natoms;
261 int min1,min2,max1,max2;
262 atom_id oindex[2];
263 rvec *x0;
264 matrix box;
265 t_trxframe frout;
266 bool bFirst;
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;
279 if (bMat) {
280 if (ng == 1) {
281 snew(leg,1);
282 sprintf(buf,"Internal in %s",grpn[0]);
283 leg[0]=strdup(buf);
284 xvgr_legend(dist,0,leg,oenv);
285 if (num) xvgr_legend(num,0,leg,oenv);
287 else {
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]);
292 leg[j]=strdup(buf);
295 xvgr_legend(dist,j,leg,oenv);
296 if (num) xvgr_legend(num,j,leg,oenv);
299 else {
300 snew(leg,ng-1);
301 for(i=0; (i<ng-1); i++) {
302 sprintf(buf,"%s-%s",grpn[0],grpn[i+1]);
303 leg[i]=strdup(buf);
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);
314 if (bPrintResName)
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");
322 j=0;
323 if (nres) {
324 snew(mindres, ng-1);
325 snew(maxdres, ng-1);
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++)
330 mindres[i-1][j]=1e6;
331 /* maxdres[*][*] is already 0 */
334 bFirst=TRUE;
335 do {
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));
344 if (bMat) {
345 if (ng == 1) {
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);
351 else {
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);
362 else {
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);
368 if (nres) {
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);
379 fprintf(dist,"\n");
380 if (num)
381 fprintf(num,"\n");
382 if ( bMin?min1:max1 != -1 )
383 if (atm)
384 fprintf(atm,"%12e %12d %12d\n",
385 output_env_conv_time(oenv,t),1+(bMin ? min1 : max1),
386 1+(bMin ? min2 : max2));
388 if (trxout>=0) {
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);
393 bFirst=FALSE;
394 /*dmin should be minimum distance for residue and group*/
395 if (bEachResEachTime)
397 fprintf(respertime,"%12e",t);
398 for(i=1; i<ng; i++)
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*/
403 mindres[i-1][j]=1e6;
404 maxdres[i-1][j]=0;
406 fprintf(respertime, "\n");
408 } while (read_next_x(oenv,status,&t,natoms,x0,box));
410 close_trj(status);
411 ffclose(dist);
412 if (num) ffclose(num);
413 if (atm) ffclose(atm);
414 if (trxout>=0) close_xtc(trxout);
416 if(nres && !bEachResEachTime) {
417 FILE *res;
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]);
427 fprintf(res, "\n");
431 sfree(x0);
434 int find_residues(t_atoms *atoms, int n, atom_id index[], atom_id **resindex)
436 int i;
437 int nres=0,resnr, presnr;
438 int *residx;
440 /* build index of first atom numbers for each residue */
441 presnr = NOTSET;
442 snew(residx, atoms->nres+1);
443 for(i=0; i<n; i++) {
444 resnr = atoms->atom[index[i]].resind;
445 if (resnr != presnr) {
446 residx[nres]=i;
447 nres++;
448 presnr=resnr;
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 */
455 residx[nres] = n;
456 *resindex = residx;
457 return nres;
460 void dump_res(FILE *out, int nres, atom_id *resindex, int n, atom_id index[])
462 int i,j;
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]);
468 fprintf(out,"\n");
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 residue 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;
501 static int ng=1;
502 static bool bEachResEachTime=FALSE,bPrintResName=FALSE;
503 t_pargs pa[] = {
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" }
525 output_env_t oenv;
526 t_topology *top=NULL;
527 int ePBC=-1;
528 char title[256];
529 real t;
530 rvec *x;
531 matrix box;
532 bool bTop=FALSE;
534 FILE *atm;
535 int i,j,nres=0;
536 const char *trxfnm,*tpsfnm,*ndxfnm,*distfnm,*numfnm,*atmfnm,*oxfnm,*resfnm;
537 char **grpname;
538 int *gnx;
539 atom_id **index, *residues=NULL;
540 t_filenm fnm[] = {
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);
567 } else {
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");
574 if (bPI) {
575 ng = 1;
576 fprintf(stderr,"Choose a group for distance calculation\n");
578 else if (!bMat)
579 ng++;
581 snew(gnx,ng);
582 snew(index,ng);
583 snew(grpname,ng);
585 if (tpsfnm || resfnm || !ndxfnm) {
586 snew(top,1);
587 bTop = read_tps_conf(tpsfnm,title,top,&ePBC,&x,NULL,box,FALSE);
588 if (bPI && !bTop)
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)) {
594 ng = gnx[0];
595 printf("Special case: making distance matrix between all atoms in group %s\n",
596 grpname[0]);
597 srenew(gnx,ng);
598 srenew(index,ng);
599 srenew(grpname,ng);
600 for(i=1; (i<ng); i++) {
601 gnx[i] = 1;
602 grpname[i] = grpname[0];
603 snew(index[i],1);
604 index[i][0] = index[0][i];
606 gnx[0] = 1;
609 if (resfnm) {
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]);
615 if (bPI) {
616 periodic_mindist_plot(trxfnm,distfnm,top,ePBC,gnx[0],index[0],bSplit,oenv);
617 } else {
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");
625 if (!bPI)
626 do_view(oenv,numfnm,"-nxy");
628 thanx(stderr);
630 return 0;