?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
changed reading hint
[gromacs/adressmacs.git] / src / kernel / gmxcheck.c
blob03cd15c20e7e9b4ec583bf7efbef599545a38725
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
11 *
12 * Copyright (c) 1991-1999
13 * BIOSON Research Institute, Dept. of Biophysical Chemistry
14 * University of Groningen, The Netherlands
15 *
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)
20 *
21 * Also check out our WWW page:
22 * http://md.chem.rug.nl/~gmx
23 * or e-mail to:
24 * gromacs@chem.rug.nl
25 *
26 * And Hey:
27 * GRowing Old MAkes el Chrono Sweat
28 */
29 static char *SRCID_gmxcheck_c = "$Id$";
30
31 #include <stdio.h>
32 #include <string.h>
33 #include <ctype.h>
34 #include "main.h"
35 #include "macros.h"
36 #include "math.h"
37 #include "futil.h"
38 #include "statutil.h"
39 #include "copyrite.h"
40 #include "sysstuff.h"
41 #include "txtdump.h"
42 #include "fatal.h"
43 #include "gmxfio.h"
44 #include "trnio.h"
45 #include "xtcio.h"
46 #include "tpbcmp.h"
47 #include "atomprop.h"
48 #include "vec.h"
49 #include "pbc.h"
50 #include "physics.h"
51 #include "smalloc.h"
52 #include "confio.h"
53 #include "enxio.h"
54
55 void chk_trj(char *fn)
56 {
57 t_trnheader sh,count;
58 int idum,j=-1,new_natoms,natoms,step;
59 real rdum,t,tt,t0,old_t1,old_t2,prec;
60 bool bShowTimestep=TRUE,bOK,newline=FALSE;
61 rvec *x;
62 matrix box;
63 size_t fpos;
64 int xd;
65 int status,ftp;
66
67 ftp = fn2ftp(fn);
68 new_natoms = -1;
69 natoms = -1;
70 t = 0;
71 t0 = NOTSET;
72 step = 0;
73 count.box_size=0;
74 count.vir_size=0;
75 count.pres_size=0;
76 count.x_size=0;
77 count.v_size=0;
78 count.f_size=0;
79
80 printf("Checking file %s\n",fn);
81 switch (ftp) {
82 case efTRJ:
83 case efTRR:
84 status = open_trn(fn,"r");
85 j = 0;
86 t = -1;
87 old_t2 = -2.0;
88 old_t1 = -1.0;
89 fpos = 0;
90 while (fread_trnheader(status,&sh,&bOK)) {
91 if (j == 0)
92 fprintf(stderr,"\n# Atoms %d\n\n",sh.natoms);
93 newline=TRUE;
94 if ((natoms > 0) && (new_natoms != natoms)) {
95 fprintf(stderr,"\nNumber of atoms at t=%g don't match (%d, %d)\n",
96 old_t1,natoms,new_natoms);
97 newline=FALSE;
98 }
99 if (j>=2) {
100 if ( fabs((sh.t-old_t1)-(old_t1-old_t2)) >
101 0.1*(fabs(sh.t-old_t1)+fabs(old_t1-old_t2)) ) {
102 bShowTimestep=FALSE;
103 fprintf(stderr,"%sTimesteps at t=%g don't match (%g, %g)\n",
104 newline?"\n":"",old_t1,old_t1-old_t2,sh.t-old_t1);
105 }
106 }
107 natoms=new_natoms;
108 old_t2=old_t1;
109 old_t1=sh.t;
110 if (t0 == NOTSET) t0=sh.t;
111 fprintf(stderr,"\rframe: %6d, step: %8d, t: %10.3f",j,sh.step,sh.t);
112 if (ftp == efTRJ)
113 fprintf(stderr," byte: %10lu",(unsigned long)fpos);
114 if (j == 0)
115 fprintf(stderr,"\n");
116 if (!fread_htrn(status,&sh,NULL,NULL,NULL,NULL))
117 fprintf(stderr,"\nframe %d at t=%g is incomplete\n",j,sh.t);
118 else {
119 j++;
120 t=sh.t;
121 new_natoms=sh.natoms;
122 #define INC(s,n,item) if (s.item != 0) n.item++
123 INC(sh,count,box_size);
124 INC(sh,count,vir_size);
125 INC(sh,count,pres_size);
126 INC(sh,count,x_size);
127 INC(sh,count,v_size);
128 INC(sh,count,f_size);
129 #undef INC
130 }
131 fpos = fio_ftell(status);
132 }
133 fprintf(stderr,"\n");
134 if (!bOK)
135 fprintf(stderr,"header of frame %d at t=%g is incomplete\n",j,sh.t);
136 close_trn(status);
137 break;
138 case efXTC:
139 xd = open_xtc(fn,"r");
140 t=-1;
141 if (read_first_xtc(xd,&new_natoms,&step,&tt,box,&x,&prec,&bOK)) {
142 fprintf(stderr,"\n# Atoms %d, XTC precision %g\n\n",new_natoms,prec);
143 j=0;
144 old_t2=-2.0;
145 old_t1=-1.0;
146 do {
147 t=tt;
148 if (j>=2) {
149 if ( fabs((t-old_t1)-(old_t1-old_t2)) >
150 0.1*(fabs(t-old_t1)+fabs(old_t1-old_t2)) ) {
151 bShowTimestep=FALSE;
152 fprintf(stderr,"%sTimesteps at t=%g don't match (%g, %g)\n",
153 newline?"\n":"",old_t1,old_t1-old_t2,t-old_t1);
154 }
155 }
156 old_t2=old_t1;
157 old_t1=t;
158 if (t0 == NOTSET) t0=t;
159 fprintf(stderr,"\rframe: %6d, step %8d, t: %10.3f",j,step,t);
160 if (j == 0)
161 fprintf(stderr,"\n");
162 newline=TRUE;
163 if ((natoms > 0) && (new_natoms != natoms)) {
164 fprintf(stderr,"\nNumber of atoms at t=%g don't match (%d, %d)\n",
165 old_t1,natoms,new_natoms);
166 newline=FALSE;
167 }
168 natoms=new_natoms;
169 count.x_size++;
170 count.box_size++;
171 j++;
172 } while (read_next_xtc(xd,&new_natoms,&step,&tt,box,x,&prec,&bOK));
173 close_xtc(xd);
174 }
175 else
176 fprintf(stderr,"Empty file %s\n",fn);
177 fprintf(stderr,"\n");
178 if (!bOK)
179 fprintf(stderr,"frame %d at t=%g is incomplete\n",j,tt);
180 break;
181 default:
182 fatal_error(0,"trajectory file .%s not supported\n",ftp2ext(ftp));
183 }
184
185 fprintf(stderr,"\nItem #frames");
186 if (bShowTimestep)
187 fprintf(stderr," Timestep (ps)");
188 fprintf(stderr,"\n");
189 #define PRINTITEM(label,item) fprintf(stderr,"%-10s %6d",label,count.item); if ((bShowTimestep) && (count.item > 1)) fprintf(stderr," %g\n",(t-t0)/(count.item-1)); else fprintf(stderr,"\n")
190 PRINTITEM ( "Box", box_size );
191 PRINTITEM ( "Virial", vir_size );
192 PRINTITEM ( "Pressure", pres_size );
193 PRINTITEM ( "Coords", x_size );
194 PRINTITEM ( "Velocities", v_size );
195 PRINTITEM ( "Forces", f_size );
196 }
197
198 void chk_tps(char *fn, real vdw_fac, real bon_lo, real bon_hi)
199 {
200 int natom,i,j,k;
201 char title[STRLEN];
202 t_topology top;
203 t_atoms *atoms;
204 rvec *x,*v;
205 rvec dx;
206 matrix box;
207 bool bV,bX,bB,bFirst,bOut;
208 real r2,ekin,temp1,temp2,dist2,vdwfac2,bonlo2,bonhi2;
209 real *atom_vdw;
210
211 fprintf(stderr,"Checking coordinate file %s\n",fn);
212 read_tps_conf(fn,title,&top,&x,&v,box,TRUE);
213 atoms=&top.atoms;
214 natom=atoms->nr;
215 fprintf(stderr,"%d atoms in file\n",atoms->nr);
216
217 /* check coordinates and box */
218 bV=FALSE;
219 bX=FALSE;
220 for (i=0; (i<natom) && !(bV && bX); i++)
221 for (j=0; (j<DIM) && !(bV && bX); j++) {
222 bV=bV || (v[i][j]!=0);
223 bX=bX || (x[i][j]!=0);
224 }
225 bB=FALSE;
226 for (i=0; (i<DIM) && !bB; i++)
227 for (j=0; (j<DIM) && !bB; j++)
228 bB=bB || (box[i][j]!=0);
229
230 fprintf(stderr,"coordinates %s\n",bX?"found":"absent");
231 fprintf(stderr,"box %s\n",bB?"found":"absent");
232 fprintf(stderr,"velocities %s\n",bV?"found":"absent");
233 fprintf(stderr,"\n");
234
235 /* check velocities */
236 if (bV) {
237 ekin=0.0;
238 for (i=0; (i<natom); i++)
239 for (j=0; (j<DIM); j++)
240 ekin+=0.5*atoms->atom[i].m*v[i][j]*v[i][j];
241 temp1=(2.0*ekin)/(natom*DIM*BOLTZ);
242 temp2=(2.0*ekin)/(natom*(DIM-1)*BOLTZ);
243 fprintf(stderr,"Kinetic energy: %g (kJ/mol)\n",ekin);
244 fprintf(stderr,"Assuming the number of degrees of freedom to be "
245 "Natoms * %d or Natoms * %d,\n"
246 "the velocities correspond to a temperature of the system\n"
247 "of %g K or %g K respectively.\n\n",DIM,DIM-1,temp1,temp2);
248 }
249
250 /* check coordinates */
251 if (bX) {
252 vdwfac2=sqr(vdw_fac);
253 bonlo2=sqr(bon_lo);
254 bonhi2=sqr(bon_hi);
255
256 fprintf(stderr,
257 "Checking for atoms closer than %g and not between %g and %g,\n"
258 "relative to sum of Van der Waals distance:\n",
259 vdw_fac,bon_lo,bon_hi);
260 snew(atom_vdw,natom);
261 for (i=0; (i<natom); i++)
262 atom_vdw[i]=get_vdw(*(atoms->resname[atoms->atom[i].resnr]),
263 *(atoms->atomname[i]),0.1)==0.0;
264
265 if (bB)
266 init_pbc(box,FALSE);
267
268 bFirst=TRUE;
269 for (i=0; (i<natom); i++) {
270 if (((i+1)%10)==0)
271 fprintf(stderr,"\r%5d",i+1);
272 for (j=i+1; (j<natom); j++) {
273 if (bB)
274 pbc_dx(x[i],x[j],dx);
275 else
276 rvec_sub(x[i],x[j],dx);
277 r2=iprod(dx,dx);
278 dist2=sqr(atom_vdw[i]+atom_vdw[j]);
279 if ( (r2<=dist2*bonlo2) ||
280 ( (r2>=dist2*bonhi2) && (r2<=dist2*vdwfac2) ) ) {
281 if (bFirst) {
282 fprintf(stderr,"\r%5s %4s %8s %5s %5s %4s %8s %5s %6s\n",
283 "atom#","name","residue","r_vdw",
284 "atom#","name","residue","r_vdw","distance");
285 bFirst=FALSE;
286 }
287 fprintf(stderr,
288 "\r%5d %4s %4s%4d %-5.3g %5d %4s %4s%4d %-5.3g %-6.4g\n",
289 i+1,*(atoms->atomname[i]),
290 *(atoms->resname[atoms->atom[i].resnr]),atoms->atom[i].resnr+1,
291 atom_vdw[i],
292 j+1,*(atoms->atomname[j]),
293 *(atoms->resname[atoms->atom[j].resnr]),atoms->atom[j].resnr+1,
294 atom_vdw[j],
295 sqrt(r2) );
296 }
297 }
298 }
299 if (bFirst)
300 fprintf(stderr,"\rno close atoms found\n");
301 fprintf(stderr,"\r \n");
302
303 if (bB) {
304 /* check box */
305 bFirst=TRUE;
306 k=0;
307 for (i=0; (i<natom) && (k<10); i++) {
308 bOut=FALSE;
309 for(j=0; (j<DIM) && !bOut; j++)
310 bOut = bOut || (x[i][j]<0) || (x[i][j]>box[j][j]);
311 if (bOut) {
312 k++;
313 if (bFirst) {
314 fprintf(stderr,"Atoms outside box ( ");
315 for (j=0; (j<DIM); j++)
316 fprintf(stderr,"%g ",box[j][j]);
317 fprintf(stderr,"):\n"
318 "(These may occur often and are normally not a problem)\n"
319 "%5s %4s %8s %5s %s\n",
320 "atom#","name","residue","r_vdw","coordinate");
321 bFirst=FALSE;
322 }
323 fprintf(stderr,
324 "%5d %4s %4s%4d %-5.3g",
325 i,*(atoms->atomname[i]),
326 *(atoms->resname[atoms->atom[i].resnr]),
327 atoms->atom[i].resnr,atom_vdw[i]);
328 for (j=0; (j<DIM); j++)
329 fprintf(stderr," %6.3g",x[i][j]);
330 fprintf(stderr,"\n");
331 }
332 }
333 if (k==10)
334 fprintf(stderr,"(maybe more)\n");
335 if (bFirst)
336 fprintf(stderr,"no atoms found outside box\n");
337 fprintf(stderr,"\n");
338 }
339 }
340 }
341
342 void chk_enx(char *fn)
343 {
344 int in,nre,frame,fnr,ndr;
345 char **enm=NULL;
346 t_energy *ee=NULL;
347 t_drblock dr;
348 bool bShowTStep;
349 real t,t0,old_t1,old_t2;
350
351 fprintf(stderr,"Checking energy file %s\n\n",fn);
352
353 in = open_enx(fn,"r");
354 do_enxnms(in,&nre,&enm);
355 fprintf(stderr,"%d groups in energy file",nre);
356 old_t2=-2.0;
357 old_t1=-1.0;
358 fnr=0;
359 t0=NOTSET;
360 bShowTStep=TRUE;
361 snew(ee,nre);
362 while (do_enx(in,&t,&frame,&nre,ee,&ndr,&dr)) {
363 if (fnr>=2) {
364 if ( fabs((t-old_t1)-(old_t1-old_t2)) >
365 0.1*(fabs(t-old_t1)+fabs(old_t1-old_t2)) ) {
366 bShowTStep=FALSE;
367 fprintf(stderr,"\nTimesteps at t=%g don't match (%g, %g)\n",
368 old_t1,old_t1-old_t2,t-old_t1);
369 }
370 }
371 old_t2=old_t1;
372 old_t1=t;
373 if (t0 == NOTSET) t0=t;
374 fprintf(stderr,"\rframe: %6d (index %6d), t: %10.3f",frame,fnr,t);
375 if (fnr == 0)
376 fprintf(stderr,"\n");
377 fnr++;
378 }
379 fprintf(stderr,"\n\nFound %d frames",fnr);
380 if (bShowTStep)
381 fprintf(stderr," with a timestep of %g ps",(t-t0)/(fnr-1));
382 fprintf(stderr,".\n");
383
384 }
385
386
387 int main(int argc,char *argv[])
388 {
389 static char *desc[] = {
390 "gmxcheck reads a trajectory ([TT].trj[tt], [TT].trr[tt] or ",
391 "[TT].xtc[tt]) or an energy file ([TT].ene[tt] or [TT].edr[tt])",
392 "and prints out useful information about them.[PAR]",
393 "For a coordinate file (generic structure file, e.g. [TT].gro[tt]) ",
394 "gmxcheck will check for presence of coordinates, velocities and box",
395 "in the file, for close contacts (smaller than [TT]-vdwfac[tt] and not",
396 "bonded, i.e. not between [TT]-bonlo[tt] and [TT]-bonhi[tt],",
397 "all relative to the",
398 "sum of both Van der Waals radii) and atoms outside the box",
399 "(these may occur often and are no problem). If velocities are present,",
400 "an estimated temperature will be calculated from them.[PAR]",
401 "The program will compare run input ([TT].tpr[tt], [TT].tpb[tt] or",
402 "[TT].tpa[tt]) files",
403 "when both [TT]-s1[tt] and [TT]-s2[tt] are supplied."
404 };
405 t_filenm fnm[] = {
406 { efTRX, "-f", NULL, ffOPTRD },
407 { efTPX, "-s1", "top1", ffOPTRD },
408 { efTPX, "-s2", "top2", ffOPTRD },
409 { efTPS, "-c", NULL, ffOPTRD },
410 { efENX, "-e", NULL, ffOPTRD },
411 { efENX, "-e1", "ener1", ffOPTRD },
412 { efENX, "-e2", "ener2", ffOPTRD }
413 };
414 #define NFILE asize(fnm)
415 char *fn1=NULL,*fn2=NULL;
416
417 static real vdw_fac=0.8;
418 static real bon_lo=0.4;
419 static real bon_hi=0.7;
420 static real ftol=0;
421 static t_pargs pa[] = {
422 { "-vdwfac", FALSE, etREAL, {&vdw_fac},
423 "Fraction of sum of VdW radii used as warning cutoff" },
424 { "-bonlo", FALSE, etREAL, {&bon_lo},
425 "Min. fract. of sum of VdW radii for bonded atoms" },
426 { "-bonhi", FALSE, etREAL, {&bon_hi},
427 "Max. fract. of sum of VdW radii for bonded atoms" },
428 { "-tol", FALSE, etREAL, {&ftol},
429 "Tolerance for comparing energy terms between different energy files" }
430 };
431
432 CopyRight(stdout,argv[0]);
433 parse_common_args(&argc,argv,0,FALSE,NFILE,fnm,asize(pa),pa,
434 asize(desc),desc,0,NULL);
435
436 if (ftp2bSet(efTRX,NFILE,fnm))
437 chk_trj(ftp2fn(efTRX,NFILE,fnm));
438
439 fn1 = opt2fn_null("-s1",NFILE,fnm);
440 fn2 = opt2fn_null("-s2",NFILE,fnm);
441 if (fn1 && fn2)
442 comp_tpx(fn1,fn2);
443 else if (fn1 || fn2)
444 fprintf(stderr,"Please give me TWO run input (.tpr/.tpa/.tpb) files!\n");
445
446 fn1 = opt2fn_null("-e1",NFILE,fnm);
447 fn2 = opt2fn_null("-e2",NFILE,fnm);
448 if (fn1 && fn2)
449 comp_enx(fn1,fn2,ftol);
450 else if (fn1 || fn2)
451 fprintf(stderr,"Please give me TWO energy (.edr/.ene) files!\n");
452
453 if (ftp2bSet(efTPS,NFILE,fnm))
454 chk_tps(ftp2fn(efTPS,NFILE,fnm), vdw_fac, bon_lo, bon_hi);
455
456 if (ftp2bSet(efENX,NFILE,fnm))
457 chk_enx(ftp2fn(efENX,NFILE,fnm));
458
459 thanx(stderr);
460
461 return 0;
462 }
got merge a long time ago