changed reading hint
[gromacs/adressmacs.git] / src / tools / g_mindist.c
blobab3519c7d2465d40e4b0656d3508044be1b63af5
1 /*
2 * $Id$
3 *
4 * This source code is part of
5 *
6 * G R O M A C S
7 *
8 * GROningen MAchine for Chemical Simulations
9 *
10 * VERSION 2.0
12 * Copyright (c) 1991-1999
13 * BIOSON Research Institute, Dept. of Biophysical Chemistry
14 * University of Groningen, The Netherlands
16 * Please refer to:
17 * GROMACS: A message-passing parallel molecular dynamics implementation
18 * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
19 * Comp. Phys. Comm. 91, 43-56 (1995)
21 * Also check out our WWW page:
22 * http://md.chem.rug.nl/~gmx
23 * or e-mail to:
24 * gromacs@chem.rug.nl
26 * And Hey:
27 * GRowing Old MAkes el Chrono Sweat
29 static char *SRCID_g_mindist_c = "$Id$";
31 #include <math.h>
32 #include <stdlib.h>
33 #include "sysstuff.h"
34 #include "string.h"
35 #include "typedefs.h"
36 #include "smalloc.h"
37 #include "macros.h"
38 #include "vec.h"
39 #include "xvgr.h"
40 #include "pbc.h"
41 #include "copyrite.h"
42 #include "futil.h"
43 #include "statutil.h"
44 #include "rdgroup.h"
46 static void calc_mindist(real MinDist,
47 matrix box, rvec x[],
48 int nx1,int nx2,
49 atom_id index1[], atom_id index2[],
50 real *md, int *nd,
51 int *ixmin, int *jxmin)
53 int i,j,j0=0,j1;
54 int ix,jx;
55 atom_id *index3;
56 rvec dx;
57 int numd=0;
58 real r2,rmin2,r_2;
60 *ixmin = -1;
61 *jxmin = -1;
63 /* real hb2; */
64 /* hb2=(sqr(box[XX][XX])+sqr(box[YY][YY])+sqr(box[ZZ][ZZ]))*0.5; */
66 r_2=sqr(MinDist);
68 /* Must init pbc every step because of pressure coupling */
69 init_pbc(box,FALSE);
70 if (index2) {
71 j0=0;
72 j1=nx2;
73 index3=index2;
75 else {
76 j1=nx1;
77 index3=index1;
80 rmin2=1000.0;
82 for(i=0; (i < nx1); i++) {
83 ix=index1[i];
84 if (!index2)
85 j0=i+1;
86 for(j=j0; (j < j1); j++) {
87 jx=index3[j];
88 if (ix != jx) {
89 pbc_dx(x[ix],x[jx],dx);
90 r2=iprod(dx,dx);
91 if (r2 < rmin2)
92 rmin2=r2;
93 if (r2 < r_2) {
94 numd++;
95 *ixmin=ix;
96 *jxmin=jx;
101 *md=sqrt(rmin2);
102 *nd=numd;
105 void mindist_plot(char *fn,FILE *atm,real mind,
106 char *dfile,char *nfile,bool bMat,
107 int ng,atom_id *index[], int gnx[], char *grpn[])
109 FILE *dist,*num;
110 char buf[256];
111 char **leg;
112 real t,md;
113 int nd,status;
114 int i,j,k,natoms;
115 int min1,min2;
116 rvec *x0;
117 matrix box;
119 if ((natoms=read_first_x(&status,fn,&t,&x0,box))==0)
120 fatal_error(0,"Could not read coordinates from statusfile\n");
122 sprintf(buf,"Number of Contacts < %g nm",mind);
123 dist=xvgropen(dfile,"Minimum Distance","Time (ps)","Distance (nm)");
124 num=xvgropen(nfile,buf,"Time (ps)","Number");
126 if (bMat) {
127 snew(leg,(ng*(ng-1))/2);
128 for(i=j=0; (i<ng-1); i++) {
129 for(k=i+1; (k<ng); k++,j++) {
130 sprintf(buf,"%s-%s",grpn[i],grpn[k]);
131 leg[j]=strdup(buf);
134 xvgr_legend(dist,j,leg);
135 xvgr_legend(num,j,leg);
137 else {
138 snew(leg,ng-1);
139 for(i=0; (i<ng-1); i++) {
140 sprintf(buf,"%s-%s",grpn[0],grpn[i+1]);
141 leg[i]=strdup(buf);
143 xvgr_legend(dist,ng-1,leg);
144 xvgr_legend(num,ng-1,leg);
146 j=0;
147 do {
148 fprintf(dist,"%12g",t);
149 fprintf(num,"%12g",t);
151 if (bMat) {
152 for(i=0; (i<ng-1); i++) {
153 for(k=i+1; (k<ng); k++) {
154 calc_mindist(mind,box,x0,gnx[i],gnx[k],
155 index[i],index[k],&md,&nd,
156 &min1,&min2);
157 fprintf(dist," %12g",md);
158 fprintf(num," %8d",nd);
162 else {
163 for(i=1; (i<ng); i++) {
164 calc_mindist(mind,box,x0,gnx[0],gnx[i],
165 index[0],index[i],&md,&nd,
166 &min1,&min2);
167 fprintf(dist," %12g",md);
168 fprintf(num," %8d",nd);
171 fprintf(dist,"\n");
172 fprintf(num,"\n");
173 if (min1 != -1)
174 fprintf(atm,"%12g %12d %12d\n",t,min1,min2);
175 } while (read_next_x(status,&t,natoms,x0,box));
177 close_trj(status);
178 fclose(dist);
179 fclose(num);
181 sfree(x0);
184 int main(int argc,char *argv[])
186 static char *desc[] = {
187 "g_mindist computes the distance between one group and a number of",
188 "other groups.",
189 "Both the smallest distance and the number of contacts within a given",
190 "distance are plotted to two separate output files"
193 static bool bMat=FALSE;
194 static real mindist=0.6;
195 t_pargs pa[] = {
196 { "-matrix", FALSE, etBOOL, {&bMat},
197 "Calculate half a matrix of group-group distances" },
198 { "-d", FALSE, etREAL, {&mindist},
199 "Distance for contacts" }
201 FILE *atm;
202 int ng;
203 char **grpname;
204 int *gnx;
205 atom_id **index;
206 t_filenm fnm[] = {
207 { efTRX, "-f", NULL, ffREAD },
208 { efNDX, NULL, NULL, ffREAD },
209 { efXVG, "-od","mindist",ffWRITE },
210 { efXVG, "-on","numcont", ffWRITE },
211 { efOUT, "-o","atm-pair", ffWRITE }
213 #define NFILE asize(fnm)
215 CopyRight(stderr,argv[0]);
216 parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME,TRUE,
217 NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL);
219 if (bMat)
220 fprintf(stderr,"You can compute all distances between a number of groups\n"
221 "How many groups do you want (>= 2) ?\n");
222 else
223 fprintf(stderr,"You can compute the distances between a first group\n"
224 "and a number of other groups.\n"
225 "How many other groups do you want (>= 1) ?\n");
226 ng=0;
227 do {
228 scanf("%d",&ng);
229 if (!bMat)
230 ng++;
231 } while (ng < 2);
232 snew(gnx,ng);
233 snew(index,ng);
234 snew(grpname,ng);
236 rd_index(ftp2fn(efNDX,NFILE,fnm),ng,gnx,index,grpname);
238 atm=ftp2FILE(efOUT,NFILE,fnm,"w");
239 mindist_plot(ftp2fn(efTRX,NFILE,fnm),atm,mindist,
240 opt2fn("-od",NFILE,fnm),opt2fn("-on",NFILE,fnm),
241 bMat,ng,index,gnx,grpname);
242 fclose(atm);
244 xvgr_file(opt2fn("-od",NFILE,fnm),"-nxy");
245 xvgr_file(opt2fn("-on",NFILE,fnm),"-nxy");
247 thanx(stdout);
249 return 0;