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
50 #include "gmx_fatal.h"
56 static void add_contact_time(int **ccount,int *ccount_nalloc,int t)
60 if (t+2 >= *ccount_nalloc) {
62 for(i=*ccount_nalloc; i<t+2; i++)
69 int gmx_dist(int argc,char *argv[])
71 const char *desc[] = {
72 "[TT]g_dist[tt] can calculate the distance between the centers of mass of two",
73 "groups of atoms as a function of time. The total distance and its",
74 "[IT]x[it]-, [IT]y[it]-, and [IT]z[it]-components are plotted.[PAR]",
75 "Or when [TT]-dist[tt] is set, print all the atoms in group 2 that are",
76 "closer than a certain distance to the center of mass of group 1.[PAR]",
77 "With options [TT]-lt[tt] and [TT]-dist[tt] the number of contacts",
78 "of all atoms in group 2 that are closer than a certain distance",
79 "to the center of mass of group 1 are plotted as a function of the time",
80 "that the contact was continously present.[PAR]",
81 "Other programs that calculate distances are [TT]g_mindist[tt]",
88 rvec *x=NULL,*v=NULL,dx;
93 int g,d,i,j,res,teller=0;
96 int ngrps; /* the number of index groups */
97 atom_id **index,max; /* the index for the atom numbers */
98 int *isize; /* the size of each group */
99 char **grpname; /* the name of each group */
102 FILE *fp=NULL,*fplt=NULL;
103 gmx_bool bCutoff,bPrintDist,bLifeTime;
105 int *contact_time=NULL,*ccount=NULL,ccount_nalloc=0,sum;
108 gmx_rmpbc_t gpbc=NULL;
110 const char *leg[4] = { "|d|","d\\sx\\N","d\\sy\\N","d\\sz\\N" };
114 static t_pargs pa[] = {
115 { "-dist", FALSE, etREAL, {&cut},
116 "Print all atoms in group 2 closer than dist to the center of mass of group 1" }
118 #define NPA asize(pa)
121 { efTRX, "-f", NULL, ffREAD },
122 { efTPX, NULL, NULL, ffREAD },
123 { efNDX, NULL, NULL, ffOPTRD },
124 { efXVG, NULL, "dist", ffOPTWR },
125 { efXVG, "-lt", "lifetime", ffOPTWR },
127 #define NFILE asize(fnm)
130 CopyRight(stderr,argv[0]);
132 parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_BE_NICE,
133 NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL,&oenv);
135 bCutoff = opt2parg_bSet("-dist",NPA,pa);
137 bLifeTime = opt2bSet("-lt",NFILE,fnm);
138 bPrintDist = (bCutoff && !bLifeTime);
140 top=read_top(ftp2fn(efTPX,NFILE,fnm),&ePBC);
142 /* read index files */
148 get_index(&top->atoms,ftp2fn(efNDX,NFILE,fnm),ngrps,isize,index,grpname);
153 for(g=0;(g<ngrps);g++) {
155 for(i=0;(i<isize[g]);i++) {
158 if (index[g][i] >= top->atoms.nr)
159 gmx_fatal(FARGS,"Atom number %d, item %d of group %d, is larger than number of atoms in the topolgy (%d)\n",index[g][i]+1,i+1,g+1,top->atoms.nr+1);
160 mass[g]+=top->atoms.atom[index[g][i]].m;
164 natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
168 gmx_fatal(FARGS,"Atom number %d in an index group is larger than number of atoms in the trajectory (%d)\n",(int)max+1,natoms);
171 /* open output file */
172 fp = xvgropen(ftp2fn(efXVG,NFILE,fnm),
173 "Distance","Time (ps)","Distance (nm)",oenv);
174 xvgr_legend(fp,4,leg,oenv);
178 snew(contact_time,isize[1]);
180 if (ePBC != epbcNONE)
185 gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box);
187 /* initialisation for correct distance calculations */
189 set_pbc(pbc,ePBC,box);
190 /* make molecules whole again */
191 gmx_rmpbc(gpbc,natoms,box,x);
193 /* calculate center of masses */
194 for(g=0;(g<ngrps);g++) {
196 copy_rvec(x[index[g][0]],com[g]);
198 for(d=0;(d<DIM);d++) {
200 for(i=0;(i<isize[g]);i++) {
201 com[g][d] += x[index[g][i]][d] * top->atoms.atom[index[g][i]].m;
203 com[g][d] /= mass[g];
209 /* write to output */
210 fprintf(fp,"%12.7f ",t);
211 for(g=0;(g<ngrps/2);g++) {
213 pbc_dx(pbc,com[2*g],com[2*g+1],dx);
215 rvec_sub(com[2*g],com[2*g+1],dx);
217 fprintf(fp,"%12.7f %12.7f %12.7f %12.7f",
218 norm(dx),dx[XX],dx[YY],dx[ZZ]);
222 for(i=0;(i<isize[1]);i++) {
225 pbc_dx(pbc,x[j],com[0],dx);
227 rvec_sub(x[j],com[0],dx);
232 res=top->atoms.atom[j].resind;
233 fprintf(stdout,"\rt: %g %d %s %d %s %g (nm)\n",
234 t,top->atoms.resinfo[res].nr,*top->atoms.resinfo[res].name,
235 j+1,*top->atoms.atomname[j],sqrt(dist2));
241 if (contact_time[i]) {
242 add_contact_time(&ccount,&ccount_nalloc,contact_time[i]-1);
251 } while (read_next_x(oenv,status,&t,natoms,x,box));
252 gmx_rmpbc_done(gpbc);
259 if (bCutoff && bLifeTime) {
260 /* Add the contacts still present in the last frame */
261 for(i=0; i<isize[1]; i++)
263 add_contact_time(&ccount,&ccount_nalloc,contact_time[i]-1);
265 sprintf(buf,"%s - %s within %g nm",
266 grpname[0],grpname[1],cut);
267 fp = xvgropen(opt2fn("-lt",NFILE,fnm),
268 buf,"Time (ps)","Number of contacts",oenv);
269 for(i=0; i<min(ccount_nalloc,teller-1); i++) {
270 /* Account for all subintervals of longer intervals */
272 for(j=i; j<ccount_nalloc; j++)
273 sum += (j-i+1)*ccount[j];
275 fprintf(fp,"%10.3f %10.3f\n",i*(t-t0)/(teller-1),sum/(double)(teller-i));