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
;
118 gmx_rmpbc_t gpbc
=NULL
;
120 natoms
=read_first_x(oenv
,&status
,trxfn
,&t
,&x
,box
);
122 check_index(NULL
,n
,index
,NULL
,natoms
);
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
);
134 gpbc
= gmx_rmpbc_init(&top
->idef
,ePBC
,natoms
,box
);
139 gmx_rmpbc(gpbc
,box
,x
,x
);
141 periodic_dist(box
,x
,n
,index
,&rmin
,&rmax
,ind_min
);
145 ind_mini
= ind_min
[0];
146 ind_minj
= ind_min
[1];
148 if ( bSplit
&& !bFirst
&& abs(t
/output_env_get_time_factor(oenv
))<1e-5 )
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]));
153 } while(read_next_x(oenv
,status
,&t
,natoms
,x
,box
));
156 gmx_rmpbc_done(gpbc
);
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);
167 static void calc_dist(real rcut
, bool bPBC
, int ePBC
, matrix box
, rvec x
[],
168 int nx1
,int nx2
, atom_id index1
[], atom_id index2
[],
170 real
*rmin
, real
*rmax
, int *nmin
, int *nmax
,
171 int *ixmin
, int *jxmin
, int *ixmax
, int *jxmax
)
177 real r2
,rmin2
,rmax2
,rcut2
;
190 /* Must init pbc every step because of pressure coupling */
192 set_pbc(&pbc
,ePBC
,box
);
205 for(j
=0; (j
< j1
); j
++) {
207 if (index2
== NULL
) {
212 for(i
=i0
; (i
< nx1
); i
++) {
216 pbc_dx(&pbc
,x
[ix
],x
[jx
],dx
);
218 rvec_sub(x
[ix
],x
[jx
],dx
);
232 } else if (r2
> rcut2
) {
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
,bool bMat
,t_atoms
*atoms
,
256 int ng
,atom_id
*index
[],int gnx
[],char *grpn
[],bool bSplit
,
257 bool bMin
, int nres
, atom_id
*residue
,bool bPBC
,int ePBC
,
258 bool bGroup
,bool bEachResEachTime
, bool bPrintResName
,
259 const output_env_t oenv
)
261 FILE *atm
,*dist
,*num
;
265 real t
,dmin
,dmax
,**mindres
=NULL
,**maxdres
=NULL
;
269 int min1
,min2
,max1
,max2
;
275 FILE *respertime
=NULL
;
277 if ((natoms
=read_first_x(oenv
,&status
,fn
,&t
,&x0
,box
))==0)
278 gmx_fatal(FARGS
,"Could not read coordinates from statusfile\n");
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
;
290 sprintf(buf
,"Internal in %s",grpn
[0]);
292 xvgr_legend(dist
,0,leg
,oenv
);
293 if (num
) xvgr_legend(num
,0,leg
,oenv
);
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
]);
303 xvgr_legend(dist
,j
,leg
,oenv
);
304 if (num
) xvgr_legend(num
,j
,leg
,oenv
);
309 for(i
=0; (i
<ng
-1); i
++) {
310 sprintf(buf
,"%s-%s",grpn
[0],grpn
[i
+1]);
313 xvgr_legend(dist
,ng
-1,leg
,oenv
);
314 if (num
) xvgr_legend(num
,ng
-1,leg
,oenv
);
317 if (bEachResEachTime
)
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,leg
,oenv
);
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");
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
++)
339 /* maxdres[*][*] is already 0 */
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");
349 fprintf(dist
,"%12e",output_env_conv_time(oenv
,t
));
350 if (num
) fprintf(num
,"%12e",output_env_conv_time(oenv
,t
));
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
);
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
);
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
);
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
,&min1
,&min2
,&max1
,&max2
);
381 mindres
[i
-1][j
] = min(mindres
[i
-1][j
],dmin
);
382 maxdres
[i
-1][j
] = max(maxdres
[i
-1][j
],dmax
);
390 if ( bMin
?min1
:max1
!= -1 )
392 fprintf(atm
,"%12e %12d %12d\n",
393 output_env_conv_time(oenv
,t
),1+(bMin
? min1
: max1
),
394 1+(bMin
? min2
: max2
));
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
);
402 /*dmin should be minimum distance for residue and group*/
403 if (bEachResEachTime
)
405 fprintf(respertime
,"%12e",t
);
407 for(j
=0; j
<nres
; j
++)
409 fprintf(respertime
, " %7g", bMin
? mindres
[i
-1][j
] : maxdres
[i
-1][j
]);
410 /*reset distances for next time point*/
414 fprintf(respertime
, "\n");
416 } while (read_next_x(oenv
,status
,&t
,natoms
,x0
,box
));
420 if (num
) ffclose(num
);
421 if (atm
) ffclose(atm
);
422 if (trxout
) close_trx(trxout
);
424 if(nres
&& !bEachResEachTime
) {
427 sprintf(buf
,"%simum Distance",bMin
? "Min" : "Max");
428 res
=xvgropen(rfile
,buf
,"Residue (#)","Distance (nm)",oenv
);
429 xvgr_legend(res
,ng
-1,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
]);
442 int find_residues(t_atoms
*atoms
, int n
, atom_id index
[], atom_id
**resindex
)
445 int nres
=0,resnr
, presnr
;
448 /* build index of first atom numbers for each residue */
450 snew(residx
, atoms
->nres
+1);
452 resnr
= atoms
->atom
[index
[i
]].resind
;
453 if (resnr
!= presnr
) {
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 */
468 void dump_res(FILE *out
, int nres
, atom_id
*resindex
, int n
, atom_id index
[])
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
]);
480 int gmx_mindist(int argc
,char *argv
[])
482 const char *desc
[] = {
483 "g_mindist 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 an other 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]."
502 const char *bugs
[] = {
503 "The [TT]-pi[tt] option is very slow."
506 static bool bMat
=FALSE
,bPI
=FALSE
,bSplit
=FALSE
,bMax
=FALSE
,bPBC
=TRUE
;
507 static bool bGroup
=FALSE
;
508 static real rcutoff
=0.6;
510 static bool bEachResEachTime
=FALSE
,bPrintResName
=FALSE
;
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" }
534 t_topology
*top
=NULL
;
544 const char *trxfnm
,*tpsfnm
,*ndxfnm
,*distfnm
,*numfnm
,*atmfnm
,*oxfnm
,*resfnm
;
547 atom_id
**index
, *residues
=NULL
;
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
}
558 #define NFILE asize(fnm)
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
);
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
);
576 tpsfnm
= ftp2fn_null(efTPS
,NFILE
,fnm
);
579 if (!tpsfnm
&& !ndxfnm
)
580 gmx_fatal(FARGS
,"You have to specify either the index file or a tpr file");
584 fprintf(stderr
,"Choose a group for distance calculation\n");
593 if (tpsfnm
|| resfnm
|| !ndxfnm
) {
595 bTop
= read_tps_conf(tpsfnm
,title
,top
,&ePBC
,&x
,NULL
,box
,FALSE
);
597 printf("\nWARNING: Without a run input file a trajectory with broken molecules will not give the correct periodic image distance\n\n");
599 get_index(top
? &(top
->atoms
) : NULL
,ndxfnm
,ng
,gnx
,index
,grpname
);
601 if (bMat
&& (ng
== 1)) {
603 printf("Special case: making distance matrix between all atoms in group %s\n",
608 for(i
=1; (i
<ng
); i
++) {
610 grpname
[i
] = grpname
[0];
612 index
[i
][0] = index
[0][i
];
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]);
624 periodic_mindist_plot(trxfnm
,distfnm
,top
,ePBC
,gnx
[0],index
[0],bSplit
,oenv
);
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
);
632 do_view(oenv
,distfnm
,"-nxy");
634 do_view(oenv
,numfnm
,"-nxy");