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"
63 int *res_ndx(t_atoms
*atoms
)
71 r0
=atoms
->atom
[0].resind
;
72 for(i
=0; (i
<atoms
->nr
); i
++)
73 rndx
[i
]=atoms
->atom
[i
].resind
-r0
;
78 int *res_natm(t_atoms
*atoms
)
85 snew(natm
,atoms
->nres
);
86 r0
=atoms
->atom
[0].resind
;
88 for(i
=0; (i
<atoms
->nres
); i
++) {
89 while ((atoms
->atom
[j
].resind
)-r0
== i
) {
98 static void calc_mat(int nres
, int natoms
, int rndx
[],
99 rvec x
[], atom_id
*index
,
100 real trunc
, real
**mdmat
,int **nmat
,int ePBC
,matrix box
)
107 set_pbc(&pbc
,ePBC
,box
);
109 for(resi
=0; (resi
<nres
); resi
++)
110 for(resj
=0; (resj
<nres
); resj
++)
111 mdmat
[resi
][resj
]=FARAWAY
;
112 for(i
=0; (i
<natoms
); i
++) {
114 for(j
=i
+1; (j
<natoms
); j
++) {
116 pbc_dx(&pbc
,x
[index
[i
]],x
[index
[j
]],ddx
);
122 mdmat
[resi
][resj
]=min(r2
,mdmat
[resi
][resj
]);
126 for(resi
=0; (resi
<nres
); resi
++) {
128 for(resj
=resi
+1; (resj
<nres
); resj
++) {
129 r
=sqrt(mdmat
[resi
][resj
]);
136 static void tot_nmat(int nres
, int natoms
, int nframes
, int **nmat
,
137 int *tot_n
, real
*mean_n
)
141 for (i
=0; (i
<nres
); i
++) {
142 for (j
=0; (j
<natoms
); j
++)
143 if (nmat
[i
][j
] != 0) {
145 mean_n
[i
]+=nmat
[i
][j
];
151 int gmx_mdmat(int argc
,char *argv
[])
153 const char *desc
[] = {
154 "[TT]g_mdmat[tt] makes distance matrices consisting of the smallest distance",
155 "between residue pairs. With [TT]-frames[tt], these distance matrices can be",
156 "stored in order to see differences in tertiary structure as a",
157 "function of time. If you choose your options unwisely, this may generate",
158 "a large output file. By default, only an averaged matrix over the whole",
159 "trajectory is output.",
160 "Also a count of the number of different atomic contacts between",
161 "residues over the whole trajectory can be made.",
162 "The output can be processed with [TT]xpm2ps[tt] to make a PostScript (tm) plot."
164 static real truncate
=1.5;
165 static gmx_bool bAtom
=FALSE
;
166 static int nlevels
=40;
168 { "-t", FALSE
, etREAL
, {&truncate
},
170 { "-nlevels", FALSE
, etINT
, {&nlevels
},
171 "Discretize distance in this number of levels" }
174 { efTRX
, "-f", NULL
, ffREAD
},
175 { efTPS
, NULL
, NULL
, ffREAD
},
176 { efNDX
, NULL
, NULL
, ffOPTRD
},
177 { efXPM
, "-mean", "dm", ffWRITE
},
178 { efXPM
, "-frames", "dmf", ffOPTWR
},
179 { efXVG
, "-no", "num",ffOPTWR
},
181 #define NFILE asize(fnm)
190 int *rndx
,*natm
,prevres
,newres
;
192 int i
,j
,nres
,natoms
,nframes
,it
,trxnat
;
195 gmx_bool bCalcN
,bFrames
;
197 char title
[256],label
[234];
200 real
**mdmat
,*resnr
,**totmdmat
;
201 int **nmat
,**totnmat
;
206 gmx_rmpbc_t gpbc
=NULL
;
208 CopyRight(stderr
,argv
[0]);
210 parse_common_args(&argc
,argv
,PCA_CAN_TIME
| PCA_BE_NICE
,NFILE
,fnm
,
211 asize(pa
),pa
,asize(desc
),desc
,0,NULL
,&oenv
);
213 fprintf(stderr
,"Will truncate at %f nm\n",truncate
);
214 bCalcN
= opt2bSet("-no",NFILE
,fnm
);
215 bFrames
= opt2bSet("-frames",NFILE
,fnm
);
217 fprintf(stderr
,"Will calculate number of different contacts\n");
219 read_tps_conf(ftp2fn(efTPS
,NFILE
,fnm
),title
,&top
,&ePBC
,&x
,NULL
,box
,FALSE
);
221 fprintf(stderr
,"Select group for analysis\n");
222 get_index(&top
.atoms
,ftp2fn_null(efNDX
,NFILE
,fnm
),1,&isize
,&index
,&grpname
);
225 snew(useatoms
.atom
,natoms
);
226 snew(useatoms
.atomname
,natoms
);
229 snew(useatoms
.resinfo
,natoms
);
231 prevres
= top
.atoms
.atom
[index
[0]].resind
;
233 for(i
=0;(i
<isize
);i
++) {
235 useatoms
.atomname
[i
]=top
.atoms
.atomname
[ii
];
236 if (top
.atoms
.atom
[ii
].resind
!= prevres
) {
237 prevres
= top
.atoms
.atom
[ii
].resind
;
239 useatoms
.resinfo
[i
] = top
.atoms
.resinfo
[prevres
];
241 fprintf(debug
,"New residue: atom %5s %5s %6d, index entry %5d, newres %5d\n",
242 *(top
.atoms
.resinfo
[top
.atoms
.atom
[ii
].resind
].name
),
243 *(top
.atoms
.atomname
[ii
]),
247 useatoms
.atom
[i
].resind
= newres
;
249 useatoms
.nres
= newres
+1;
252 rndx
=res_ndx(&(useatoms
));
253 natm
=res_natm(&(useatoms
));
255 fprintf(stderr
,"There are %d residues with %d atoms\n",nres
,natoms
);
263 for(i
=0; (i
<nres
); i
++) {
265 snew(nmat
[i
],natoms
);
266 snew(totnmat
[i
],natoms
);
270 for(i
=0; (i
<nres
); i
++)
271 snew(totmdmat
[i
],nres
);
273 trxnat
=read_first_x(oenv
,&status
,ftp2fn(efTRX
,NFILE
,fnm
),&t
,&x
,box
);
277 rlo
.r
=1.0, rlo
.g
=1.0, rlo
.b
=1.0;
278 rhi
.r
=0.0, rhi
.g
=0.0, rhi
.b
=0.0;
280 gpbc
= gmx_rmpbc_init(&top
.idef
,ePBC
,trxnat
,box
);
283 out
=opt2FILE("-frames",NFILE
,fnm
,"w");
285 gmx_rmpbc(gpbc
,trxnat
,box
,x
);
287 calc_mat(nres
,natoms
,rndx
,x
,index
,truncate
,mdmat
,nmat
,ePBC
,box
);
288 for (i
=0; (i
<nres
); i
++)
289 for (j
=0; (j
<natoms
); j
++)
292 for (i
=0; (i
<nres
); i
++)
293 for (j
=0; (j
<nres
); j
++)
294 totmdmat
[i
][j
] += mdmat
[i
][j
];
296 sprintf(label
,"t=%.0f ps",t
);
297 write_xpm(out
,0,label
,"Distance (nm)","Residue Index","Residue Index",
298 nres
,nres
,resnr
,resnr
,mdmat
,0,truncate
,rlo
,rhi
,&nlevels
);
300 } while (read_next_x(oenv
,status
,&t
,trxnat
,x
,box
));
301 fprintf(stderr
,"\n");
303 gmx_rmpbc_done(gpbc
);
307 fprintf(stderr
,"Processed %d frames\n",nframes
);
309 for (i
=0; (i
<nres
); i
++)
310 for (j
=0; (j
<nres
); j
++)
311 totmdmat
[i
][j
] /= nframes
;
312 write_xpm(opt2FILE("-mean",NFILE
,fnm
,"w"),0,"Mean smallest distance",
313 "Distance (nm)","Residue Index","Residue Index",
314 nres
,nres
,resnr
,resnr
,totmdmat
,0,truncate
,rlo
,rhi
,&nlevels
);
317 tot_nmat(nres
,natoms
,nframes
,totnmat
,tot_n
,mean_n
);
318 fp
=xvgropen(ftp2fn(efXVG
,NFILE
,fnm
),
319 "Increase in number of contacts","Residue","Ratio",oenv
);
320 fprintf(fp
,"@ legend on\n");
321 fprintf(fp
,"@ legend box on\n");
322 fprintf(fp
,"@ legend loctype view\n");
323 fprintf(fp
,"@ legend 0.75, 0.8\n");
324 fprintf(fp
,"@ legend string 0 \"Total/mean\"\n");
325 fprintf(fp
,"@ legend string 1 \"Total\"\n");
326 fprintf(fp
,"@ legend string 2 \"Mean\"\n");
327 fprintf(fp
,"@ legend string 3 \"# atoms\"\n");
328 fprintf(fp
,"@ legend string 4 \"Mean/# atoms\"\n");
329 fprintf(fp
,"#%3s %8s %3s %8s %3s %8s\n",
330 "res","ratio","tot","mean","natm","mean/atm");
331 for (i
=0; (i
<nres
); i
++) {
335 ratio
=tot_n
[i
]/mean_n
[i
];
336 fprintf(fp
,"%3d %8.3f %3d %8.3f %3d %8.3f\n",
337 i
+1,ratio
,tot_n
[i
],mean_n
[i
],natm
[i
],mean_n
[i
]/natm
[i
]);