?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Replaced all occurences of str(n)casecmp with gmx_str(n)casecmp.
[gromacs/rigid-bodies.git] / src / gmxlib / trxio.c
blob4bd9f29d3adb39ffe7fc5a28e62b0b7ae39d07ca
1 /*
2 *
3 * This source code is part of
4 *
5 * G R O M A C S
6 *
7 * GROningen MAchine for Chemical Simulations
8 *
9 * VERSION 3.2.0
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.
14
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.
19 *
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 forxd
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.
26 *
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.
29 *
30 * For more info, check our website at http://www.gromacs.org
31 *
32 * And Hey:
33 * GROningen Mixture of Alchemy and Childrens' Stories
34 */
35 #ifdef HAVE_CONFIG_H
36 #include <config.h>
37 #endif
38
39 #include <ctype.h>
40 #include "sysstuff.h"
41 #include "string2.h"
42 #include "smalloc.h"
43 #include "pbc.h"
44 #include "statutil.h"
45 #include "gmxfio.h"
46 #include "trnio.h"
47 #include "names.h"
48 #include "vec.h"
49 #include "futil.h"
50 #include "gmxfio.h"
51 #include "xtcio.h"
52 #include "pdbio.h"
53 #include "confio.h"
54 #include "checkpoint.h"
55 #include "wgms.h"
56 #include "vmdio.h"
57 #include <math.h>
58
59 /* defines for frame counter output */
60 #define SKIP1 10
61 #define SKIP2 100
62 #define SKIP3 1000
63
64 /* Globals for gromos-87 input */
65 typedef enum { effXYZ, effXYZBox, effG87, effG87Box, effNR } eFileFormat;
66
67 struct t_trxstatus
68 {
69 int __frame;
70 t_trxframe *xframe;
71 int nxframe;
72 t_fileio *fio;
73 eFileFormat eFF;
74 int NATOMS;
75 double DT,BOX[3];
76 bool bReadBox;
77 };
78
79 static void initcount(t_trxstatus *status)
80 {
81 status->__frame=-1;
82 }
83
84 static void status_init(t_trxstatus *status)
85 {
86 status->nxframe=0;
87 status->xframe=NULL;
88 status->fio=NULL;
89 status->__frame=-1;
90 }
91
92
93 int nframes_read(t_trxstatus *status)
94 {
95 return status->__frame;
96 }
97
98 static void printcount_(t_trxstatus *status, const output_env_t oenv,
99 const char *l,real t)
100 {
101 if ((status->__frame < 2*SKIP1 || status->__frame % SKIP1 == 0) &&
102 (status->__frame < 2*SKIP2 || status->__frame % SKIP2 == 0) &&
103 (status->__frame < 2*SKIP3 || status->__frame % SKIP3 == 0))
104 fprintf(stderr,"\r%-14s %6d time %8.3f ",l,status->__frame,
105 output_env_conv_time(oenv,t));
106 }
107
108 static void printcount(t_trxstatus *status, const output_env_t oenv,real t,
109 bool bSkip)
110 {
111 status->__frame++;
112 printcount_(status, oenv,bSkip ? "Skipping frame" : "Reading frame",t);
113 }
114
115 static void printlast(t_trxstatus *status, const output_env_t oenv,real t)
116 {
117 printcount_(status, oenv,"Last frame",t);
118 fprintf(stderr,"\n");
119 }
120
121 static void printincomp(t_trxstatus *status, t_trxframe *fr)
122 {
123 if (fr->not_ok & HEADER_NOT_OK)
124 fprintf(stderr,"WARNING: Incomplete header: nr %d time %g\n",
125 status->__frame+1,fr->time);
126 else if (fr->not_ok)
127 fprintf(stderr,"WARNING: Incomplete frame: nr %d time %g\n",
128 status->__frame+1,fr->time);
129 }
130
131 int prec2ndec(real prec)
132 {
133 if (prec <= 0)
134 gmx_fatal(FARGS,"DEATH HORROR prec (%g) <= 0 in prec2ndec",prec);
135
136 return (int)(log(prec)/log(10.0)+0.5);
137 }
138
139
140 t_fileio *trx_get_fileio(t_trxstatus *status)
141 {
142 return status->fio;
143 }
144
145
146
147 void clear_trxframe(t_trxframe *fr,bool bFirst)
148 {
149 fr->not_ok = 0;
150 fr->bTitle = FALSE;
151 fr->bStep = FALSE;
152 fr->bTime = FALSE;
153 fr->bLambda = FALSE;
154 fr->bAtoms = FALSE;
155 fr->bPrec = FALSE;
156 fr->bX = FALSE;
157 fr->bV = FALSE;
158 fr->bF = FALSE;
159 fr->bBox = FALSE;
160 if (bFirst) {
161 fr->flags = 0;
162 fr->bDouble= FALSE;
163 fr->natoms = -1;
164 fr->t0 = 0;
165 fr->tpf = 0;
166 fr->tppf = 0;
167 fr->title = NULL;
168 fr->step = 0;
169 fr->time = 0;
170 fr->lambda = 0;
171 fr->atoms = NULL;
172 fr->prec = 0;
173 fr->x = NULL;
174 fr->v = NULL;
175 fr->f = NULL;
176 clear_mat(fr->box);
177 fr->bPBC = FALSE;
178 fr->ePBC = -1;
179 }
180 }
181
182 void set_trxframe_ePBC(t_trxframe *fr,int ePBC)
183 {
184 fr->bPBC = (ePBC == -1);
185 fr->ePBC = ePBC;
186 }
187
188 int write_trxframe_indexed(t_trxstatus *status,t_trxframe *fr,int nind,
189 atom_id *ind, gmx_conect gc)
190 {
191 char title[STRLEN];
192 rvec *xout=NULL,*vout=NULL,*fout=NULL;
193 int i;
194 real prec;
195
196 if (fr->bPrec)
197 prec = fr->prec;
198 else
199 prec = 1000.0;
200
201 switch (gmx_fio_getftp(status->fio)) {
202 case efTRJ:
203 case efTRR:
204 break;
205 default:
206 if (!fr->bX)
207 gmx_fatal(FARGS,"Need coordinates to write a %s trajectory",
208 ftp2ext(gmx_fio_getftp(status->fio)));
209 break;
210 }
211
212 switch (gmx_fio_getftp(status->fio)) {
213 case efTRJ:
214 case efTRR:
215 if (fr->bV) {
216 snew(vout,nind);
217 for(i=0; i<nind; i++)
218 copy_rvec(fr->v[ind[i]],vout[i]);
219 }
220 if (fr->bF) {
221 snew(fout,nind);
222 for(i=0; i<nind; i++)
223 copy_rvec(fr->f[ind[i]],fout[i]);
224 }
225 case efXTC:
226 case efG87:
227 if (fr->bX) {
228 snew(xout,nind);
229 for(i=0; i<nind; i++)
230 copy_rvec(fr->x[ind[i]],xout[i]);
231 }
232 break;
233 default:
234 break;
235 }
236
237 switch (gmx_fio_getftp(status->fio)) {
238 case efXTC:
239 write_xtc(status->fio,nind,fr->step,fr->time,fr->box,xout,prec);
240 break;
241 case efTRJ:
242 case efTRR:
243 fwrite_trn(status->fio,nframes_read(status),
244 fr->time,fr->step,fr->box,nind,xout,vout,fout);
245 break;
246 case efGRO:
247 case efPDB:
248 case efBRK:
249 case efENT:
250 if (!fr->bAtoms)
251 gmx_fatal(FARGS,"Can not write a %s file without atom names",
252 ftp2ext(gmx_fio_getftp(status->fio)));
253 sprintf(title,"frame t= %.3f",fr->time);
254 if (gmx_fio_getftp(status->fio) == efGRO)
255 write_hconf_indexed_p(gmx_fio_getfp(status->fio),title,fr->atoms,nind,ind,
256 prec2ndec(prec),
257 fr->x,fr->bV ? fr->v : NULL,fr->box);
258 else
259 write_pdbfile_indexed(gmx_fio_getfp(status->fio),title,fr->atoms,
260 fr->x,-1,fr->box,' ',fr->step,nind,ind,gc,TRUE);
261 break;
262 case efG87:
263 write_gms(gmx_fio_getfp(status->fio),nind,xout,fr->box);
264 break;
265 case efG96:
266 write_g96_conf(gmx_fio_getfp(status->fio),fr,nind,ind);
267 break;
268 default:
269 gmx_fatal(FARGS,"Sorry, write_trxframe_indexed can not write %s",
270 ftp2ext(gmx_fio_getftp(status->fio)));
271 break;
272 }
273
274 switch (gmx_fio_getftp(status->fio)) {
275 case efTRN:
276 case efTRJ:
277 case efTRR:
278 if (vout) sfree(vout);
279 if (fout) sfree(fout);
280 case efXTC:
281 case efG87:
282 sfree(xout);
283 break;
284 default:
285 break;
286 }
287
288 return 0;
289 }
290
291 int write_trxframe(t_trxstatus *status,t_trxframe *fr,gmx_conect gc)
292 {
293 char title[STRLEN];
294 real prec;
295
296 if (fr->bPrec)
297 prec = fr->prec;
298 else
299 prec = 1000.0;
300
301 switch (gmx_fio_getftp(status->fio)) {
302 case efTRJ:
303 case efTRR:
304 break;
305 default:
306 if (!fr->bX)
307 gmx_fatal(FARGS,"Need coordinates to write a %s trajectory",
308 ftp2ext(gmx_fio_getftp(status->fio)));
309 break;
310 }
311
312 switch (gmx_fio_getftp(status->fio)) {
313 case efXTC:
314 write_xtc(status->fio,fr->natoms,fr->step,fr->time,fr->box,fr->x,prec);
315 break;
316 case efTRJ:
317 case efTRR:
318 fwrite_trn(status->fio,fr->step,fr->time,fr->lambda,fr->box,fr->natoms,
319 fr->bX ? fr->x:NULL,fr->bV ? fr->v:NULL ,fr->bF ? fr->f:NULL);
320 break;
321 case efGRO:
322 case efPDB:
323 case efBRK:
324 case efENT:
325 if (!fr->bAtoms)
326 gmx_fatal(FARGS,"Can not write a %s file without atom names",
327 ftp2ext(gmx_fio_getftp(status->fio)));
328 sprintf(title,"frame t= %.3f",fr->time);
329 if (gmx_fio_getftp(status->fio) == efGRO)
330 write_hconf_p(gmx_fio_getfp(status->fio),title,fr->atoms,
331 prec2ndec(prec),fr->x,fr->bV ? fr->v : NULL,fr->box);
332 else
333 write_pdbfile(gmx_fio_getfp(status->fio),title,
334 fr->atoms,fr->x,fr->bPBC ? fr->ePBC : -1,fr->box,
335 ' ',fr->step,gc,TRUE);
336 break;
337 case efG87:
338 write_gms(gmx_fio_getfp(status->fio),fr->natoms,fr->x,fr->box);
339 break;
340 case efG96:
341 write_g96_conf(gmx_fio_getfp(status->fio),fr,-1,NULL);
342 break;
343 default:
344 gmx_fatal(FARGS,"Sorry, write_trxframe can not write %s",
345 ftp2ext(gmx_fio_getftp(status->fio)));
346 break;
347 }
348
349 return 0;
350 }
351
352 int write_trx(t_trxstatus *status,int nind,atom_id *ind,t_atoms *atoms,
353 int step,real time,matrix box,rvec x[],rvec *v,
354 gmx_conect gc)
355 {
356 t_trxframe fr;
357
358 clear_trxframe(&fr,TRUE);
359 fr.bStep = TRUE;
360 fr.step = step;
361 fr.bTime = TRUE;
362 fr.time = time;
363 fr.bAtoms = atoms!=NULL;
364 fr.atoms = atoms;
365 fr.bX = TRUE;
366 fr.x = x;
367 fr.bV = v!=NULL;
368 fr.v = v;
369 fr.bBox = TRUE;
370 copy_mat(box,fr.box);
371
372 return write_trxframe_indexed(status,&fr,nind,ind,gc);
373 }
374
375 void close_trx(t_trxstatus *status)
376 {
377 gmx_fio_close(status->fio);
378 sfree(status);
379 }
380
381 t_trxstatus *open_trx(const char *outfile,const char *filemode)
382 {
383 t_trxstatus *stat;
384 if (filemode[0]!='w' && filemode[0]!='a' && filemode[1]!='+')
385 gmx_fatal(FARGS,"Sorry, write_trx can only write");
386
387 snew(stat,1);
388 status_init(stat);
389
390 stat->fio=gmx_fio_open(outfile,filemode);
391 return stat;
392 }
393
394 static bool gmx_next_frame(t_trxstatus *status,t_trxframe *fr)
395 {
396 t_trnheader sh;
397 bool bOK,bRet;
398
399 bRet = FALSE;
400
401 if (fread_trnheader(status->fio,&sh,&bOK)) {
402 fr->bDouble=sh.bDouble;
403 fr->natoms=sh.natoms;
404 fr->bStep=TRUE;
405 fr->step=sh.step;
406 fr->bTime=TRUE;
407 fr->time=sh.t;
408 fr->bLambda = TRUE;
409 fr->lambda = sh.lambda;
410 fr->bBox = sh.box_size>0;
411 if (fr->flags & (TRX_READ_X | TRX_NEED_X)) {
412 if (fr->x==NULL)
413 snew(fr->x,sh.natoms);
414 fr->bX = sh.x_size>0;
415 }
416 if (fr->flags & (TRX_READ_V | TRX_NEED_V)) {
417 if (fr->v==NULL)
418 snew(fr->v,sh.natoms);
419 fr->bV = sh.v_size>0;
420 }
421 if (fr->flags & (TRX_READ_F | TRX_NEED_F)) {
422 if (fr->f==NULL)
423 snew(fr->f,sh.natoms);
424 fr->bF = sh.f_size>0;
425 }
426 if (fread_htrn(status->fio,&sh,fr->box,fr->x,fr->v,fr->f))
427 bRet = TRUE;
428 else
429 fr->not_ok = DATA_NOT_OK;
430 } else
431 if (!bOK)
432 fr->not_ok = HEADER_NOT_OK;
433
434 return bRet;
435 }
436
437 static void choose_ff(FILE *fp)
438 {
439 int i,m,c;
440 int rc;
441 eFileFormat eFF;
442 t_trxstatus *stat;
443
444 printf("\n\n");
445 printf(" Select File Format\n");
446 printf("---------------------------\n");
447 printf("1. XYZ File\n");
448 printf("2. XYZ File with Box\n");
449 printf("3. Gromos-87 Ascii Trajectory\n");
450 printf("4. Gromos-87 Ascii Trajectory with Box\n");
451
452 snew(stat,1);
453 status_init(stat);
454
455 do {
456 printf("\nChoice: ");
457 fflush(stdout);
458 do
459 {
460 rc = scanf("%d",&i);
461 }
462 while (rc!=1);
463 i--;
464 } while ((i < 0) || (i >= effNR));
465 printf("\n");
466
467 stat->eFF = (eFileFormat) i;
468
469 for(m=0; (m<DIM); m++) stat->BOX[m]=0;
470
471 stat->bReadBox = (stat->eFF == effG87Box) || (stat->eFF == effXYZBox);
472
473 switch (stat->eFF) {
474 case effXYZ:
475 case effXYZBox:
476 if( 5 != fscanf(fp,"%d%lf%lf%lf%lf",&stat->NATOMS,&stat->BOX[XX],&stat->BOX[YY],&stat->BOX[ZZ],&stat->DT))
477 {
478 gmx_fatal(FARGS,"Error reading natoms/box in file");
479 }
480 break;
481 case effG87:
482 case effG87Box:
483 printf("GROMOS! OH DEAR...\n\n");
484 printf("Number of atoms ? ");
485 fflush(stdout);
486 if (1 != scanf("%d",&stat->NATOMS))
487 {
488 gmx_fatal(FARGS,"Error reading natoms in file");
489 }
490
491 printf("Time between timeframes ? ");
492 fflush(stdout);
493 if( 1 != scanf("%lf",&stat->DT))
494 {
495 gmx_fatal(FARGS,"Error reading dt from file");
496 }
497
498 if (stat->eFF == effG87) {
499 printf("Box X Y Z ? ");
500 fflush(stdout);
501 if(3 != scanf("%lf%lf%lf",&stat->BOX[XX],&stat->BOX[YY],&stat->BOX[ZZ]))
502 {
503 gmx_fatal(FARGS,"Error reading box in file");
504 }
505 }
506 do {
507 c=fgetc(fp);
508 printf("%c",c);
509 } while (c != '\n');
510 printf("\n");
511 fflush(stdout);
512 break;
513 default:
514 printf("Hellow World\n");
515 }
516 }
517
518 static bool do_read_xyz(t_trxstatus *status, FILE *fp,int natoms,
519 rvec x[],matrix box)
520 {
521 int i,m;
522 double x0;
523
524 for(i=0; (i<natoms); i++) {
525 for(m=0; (m<DIM); m++) {
526 if (fscanf(fp,"%lf",&x0) != 1) {
527 if (i || m)
528 fprintf(stderr,"error reading statusfile: x[%d][%d]\n",i,m);
529 /* else eof! */
530 return FALSE;
531 }
532 x[i][m]=x0;
533 }
534 }
535 if (status->bReadBox) {
536 for(m=0; (m<DIM); m++) {
537 if (fscanf(fp,"%lf",&x0) != 1)
538 return FALSE;
539 box[m][m]=x0;
540 }
541 }
542 return TRUE;
543 }
544
545 static bool xyz_next_x(t_trxstatus *status, FILE *fp, const output_env_t oenv,
546 real *t, int natoms, rvec x[], matrix box)
547 /* Reads until a new x can be found (return TRUE)
548 * or eof (return FALSE)
549 */
550 {
551 real pt;
552
553 pt=*t;
554 while (!bTimeSet(TBEGIN) || (*t < rTimeValue(TBEGIN))) {
555 if (!do_read_xyz(status,fp,natoms,x,box))
556 return FALSE;
557 printcount(status,oenv,*t,FALSE);
558 *t+=status->DT;
559 pt=*t;
560 }
561 if (!bTimeSet(TEND) || (*t <= rTimeValue(TEND))) {
562 if (!do_read_xyz(status,fp,natoms,x,box)) {
563 printlast(status, oenv,*t);
564 return FALSE;
565 }
566 printcount(status,oenv,*t,FALSE);
567 pt=*t;
568 *t+=status->DT;
569 return TRUE;
570 }
571 printlast(status,oenv,pt);
572 return FALSE;
573 }
574
575 static int xyz_first_x(t_trxstatus *status, FILE *fp, const output_env_t oenv,
576 real *t, rvec **x, matrix box)
577 /* Reads fp, mallocs x, and returns x and box
578 * Returns natoms when successful, FALSE otherwise
579 */
580 {
581 int m;
582
583 initcount(status);
584
585 clear_mat(box);
586 choose_ff(fp);
587
588 for(m=0; (m<DIM); m++)
589 box[m][m]=status->BOX[m];
590
591 snew(*x,status->NATOMS);
592 *t=status->DT;
593 if (!xyz_next_x(status, fp,oenv,t,status->NATOMS,*x,box))
594 return 0;
595 *t=0.0;
596
597 return status->NATOMS;
598 }
599
600 static bool pdb_next_x(t_trxstatus *status, FILE *fp,t_trxframe *fr)
601 {
602 t_atoms atoms;
603 matrix boxpdb;
604 int ePBC,model_nr,na;
605 char title[STRLEN],*time;
606 double dbl;
607
608 atoms.nr = fr->natoms;
609 atoms.atom=NULL;
610 atoms.pdbinfo=NULL;
611 /* the other pointers in atoms should not be accessed if these are NULL */
612 model_nr=NOTSET;
613 na=read_pdbfile(fp,title,&model_nr,&atoms,fr->x,&ePBC,boxpdb,TRUE,NULL);
614 set_trxframe_ePBC(fr,ePBC);
615 if (nframes_read(status)==0)
616 fprintf(stderr," '%s', %d atoms\n",title, fr->natoms);
617 fr->bPrec = TRUE;
618 fr->prec = 10000;
619 fr->bX = TRUE;
620 fr->bBox = (boxpdb[XX][XX] != 0.0);
621 if (fr->bBox) {
622 copy_mat(boxpdb,fr->box);
623 }
624
625 if (model_nr!=NOTSET) {
626 fr->bStep = TRUE;
627 fr->step = model_nr;
628 }
629 time=strstr(title," t= ");
630 if (time) {
631 fr->bTime = TRUE;
632 sscanf(time+4,"%lf",&dbl);
633 fr->time=(real)dbl;
634 } else {
635 fr->bTime = FALSE;
636 /* this is a bit dirty, but it will work: if no time is read from
637 comment line in pdb file, set time to current frame number */
638 if (fr->bStep)
639 fr->time=(real)fr->step;
640 else
641 fr->time=(real)nframes_read(status);
642 }
643 if (na == 0) {
644 return FALSE;
645 } else {
646 if (na != fr->natoms)
647 gmx_fatal(FARGS,"Number of atoms in pdb frame %d is %d instead of %d",
648 nframes_read(status),na,fr->natoms);
649 return TRUE;
650 }
651 }
652
653 static int pdb_first_x(t_trxstatus *status, FILE *fp, t_trxframe *fr)
654 {
655 initcount(status);
656
657 fprintf(stderr,"Reading frames from pdb file");
658 frewind(fp);
659 get_pdb_coordnum(fp, &fr->natoms);
660 if (fr->natoms==0)
661 gmx_fatal(FARGS,"\nNo coordinates in pdb file\n");
662 frewind(fp);
663 snew(fr->x,fr->natoms);
664 pdb_next_x(status, fp, fr);
665
666 return fr->natoms;
667 }
668
669 bool read_next_frame(const output_env_t oenv,t_trxstatus *status,t_trxframe *fr)
670 {
671 real pt;
672 int ct;
673 bool bOK,bRet,bMissingData=FALSE,bSkip=FALSE;
674 int dummy=0;
675
676 bRet = FALSE;
677 pt=fr->time;
678
679 do {
680 clear_trxframe(fr,FALSE);
681 fr->tppf = fr->tpf;
682 fr->tpf = fr->time;
683
684 switch (gmx_fio_getftp(status->fio)) {
685 case efTRJ:
686 case efTRR:
687 bRet = gmx_next_frame(status,fr);
688 break;
689 case efCPT:
690 /* Checkpoint files can not contain mulitple frames */
691 break;
692 case efG96:
693 read_g96_conf(gmx_fio_getfp(status->fio),NULL,fr);
694 bRet = (fr->natoms > 0);
695 break;
696 case efG87:
697 bRet = xyz_next_x(status, gmx_fio_getfp(status->fio),oenv,&fr->time,
698 fr->natoms, fr->x,fr->box);
699 fr->bTime = bRet;
700 fr->bX = bRet;
701 fr->bBox = bRet;
702 break;
703 case efXTC:
704 /* B. Hess 2005-4-20
705 * Sometimes is off by one frame
706 * and sometimes reports frame not present/file not seekable
707 */
708 /* DvdS 2005-05-31: this has been fixed along with the increased
709 * accuracy of the control over -b and -e options.
710 */
711 if (bTimeSet(TBEGIN) && (fr->time < rTimeValue(TBEGIN))) {
712 if (xtc_seek_time(status->fio, rTimeValue(TBEGIN),fr->natoms)) {
713 gmx_fatal(FARGS,"Specified frame doesn't exist or file not seekable");
714 }
715 initcount(status);
716 }
717 bRet = read_next_xtc(status->fio,fr->natoms,&fr->step,&fr->time,fr->box,
718 fr->x,&fr->prec,&bOK);
719 fr->bPrec = (bRet && fr->prec > 0);
720 fr->bStep = bRet;
721 fr->bTime = bRet;
722 fr->bX = bRet;
723 fr->bBox = bRet;
724 if (!bOK) {
725 /* Actually the header could also be not ok,
726 but from bOK from read_next_xtc this can't be distinguished */
727 fr->not_ok = DATA_NOT_OK;
728 }
729 break;
730 case efPDB:
731 bRet = pdb_next_x(status, gmx_fio_getfp(status->fio),fr);
732 break;
733 case efGRO:
734 bRet = gro_next_x_or_v(gmx_fio_getfp(status->fio),fr);
735 break;
736 default:
737 #ifdef GMX_DLOPEN
738 bRet = read_next_vmd_frame(dummy,fr);
739 #else
740 gmx_fatal(FARGS,"DEATH HORROR in read_next_frame ftp=%s,status=%s",
741 ftp2ext(gmx_fio_getftp(status->fio)),
742 gmx_fio_getname(status->fio));
743 #endif
744 }
745
746 if (bRet) {
747 bMissingData = ((fr->flags & TRX_NEED_X && !fr->bX) ||
748 (fr->flags & TRX_NEED_V && !fr->bV) ||
749 (fr->flags & TRX_NEED_F && !fr->bF));
750 bSkip = FALSE;
751 if (!bMissingData) {
752 ct=check_times2(fr->time,fr->t0,fr->tpf,fr->tppf,fr->bDouble);
753 if (ct == 0 || (fr->flags & TRX_DONT_SKIP && ct<0)) {
754 printcount(status, oenv,fr->time,FALSE);
755 } else if (ct > 0)
756 bRet = FALSE;
757 else {
758 printcount(status, oenv,fr->time,TRUE);
759 bSkip = TRUE;
760 }
761 }
762 }
763
764 } while (bRet && (bMissingData || bSkip));
765
766 if (!bRet) {
767 printlast(status, oenv,pt);
768 if (fr->not_ok)
769 printincomp(status, fr);
770 }
771
772 return bRet;
773 }
774
775 int read_first_frame(const output_env_t oenv,t_trxstatus **status,
776 const char *fn,t_trxframe *fr,int flags)
777 {
778 t_fileio *fio;
779 bool bFirst,bOK;
780 int dummy=0;
781
782 clear_trxframe(fr,TRUE);
783 fr->flags = flags;
784
785 bFirst = TRUE;
786
787 snew((*status), 1);
788
789 status_init( *status );
790 (*status)->nxframe=1;
791 initcount(*status);
792
793 fio = (*status)->fio =gmx_fio_open(fn,"r");
794 switch (gmx_fio_getftp(fio))
795 {
796 case efTRJ:
797 case efTRR:
798 break;
799 case efCPT:
800 read_checkpoint_trxframe(fio,fr);
801 bFirst = FALSE;
802 break;
803 case efG96:
804 /* Can not rewind a compressed file, so open it twice */
805 read_g96_conf(gmx_fio_getfp(fio),fn,fr);
806 gmx_fio_close(fio);
807 clear_trxframe(fr,FALSE);
808 if (flags & (TRX_READ_X | TRX_NEED_X))
809 snew(fr->x,fr->natoms);
810 if (flags & (TRX_READ_V | TRX_NEED_V))
811 snew(fr->v,fr->natoms);
812 fio = (*status)->fio =gmx_fio_open(fn,"r");
813 break;
814 case efG87:
815 fr->natoms=xyz_first_x(*status, gmx_fio_getfp(fio),oenv,&fr->time,
816 &fr->x,fr->box);
817 if (fr->natoms) {
818 fr->bTime = TRUE;
819 fr->bX = TRUE;
820 fr->bBox = TRUE;
821 printcount(*status,oenv,fr->time,FALSE);
822 }
823 bFirst = FALSE;
824 break;
825 case efXTC:
826 if (read_first_xtc(fio,&fr->natoms,&fr->step,&fr->time,fr->box,&fr->x,
827 &fr->prec,&bOK) == 0) {
828 if (bOK) {
829 gmx_fatal(FARGS,"No XTC!\n");
830 } else {
831 fr->not_ok = DATA_NOT_OK;
832 }
833 }
834 if (fr->not_ok) {
835 fr->natoms = 0;
836 printincomp(*status,fr);
837 } else {
838 fr->bPrec = (fr->prec > 0);
839 fr->bStep = TRUE;
840 fr->bTime = TRUE;
841 fr->bX = TRUE;
842 fr->bBox = TRUE;
843 printcount(*status,oenv,fr->time,FALSE);
844 }
845 bFirst = FALSE;
846 break;
847 case efPDB:
848 pdb_first_x(*status, gmx_fio_getfp(fio),fr);
849 if (fr->natoms)
850 printcount(*status,oenv,fr->time,FALSE);
851 bFirst = FALSE;
852 break;
853 case efGRO:
854 if (gro_first_x_or_v(gmx_fio_getfp(fio),fr))
855 printcount(*status,oenv,fr->time,FALSE);
856 bFirst = FALSE;
857 break;
858 default:
859 #ifdef GMX_DLOPEN
860 gmx_fio_fp_close(fio); /*only close the file without removing FIO entry*/
861 if (!read_first_vmd_frame(&dummy,fn,fr,flags))
862 {
863 gmx_fatal(FARGS,"Not supported in read_first_frame: %s",fn);
864 }
865 #else
866 gmx_fatal(FARGS,"Not supported in read_first_frame: %s",fn);
867 #endif
868 break;
869 }
870
871 /* Return FALSE if we read a frame that's past the set ending time. */
872 if (!bFirst && (!(fr->flags & TRX_DONT_SKIP) && check_times(fr->time) > 0)) {
873 fr->t0 = fr->time;
874 return FALSE;
875 }
876
877 if (bFirst ||
878 (!(fr->flags & TRX_DONT_SKIP) && check_times(fr->time) < 0))
879 /* Read a frame when no frame was read or the first was skipped */
880 if (!read_next_frame(oenv,*status,fr))
881 return FALSE;
882 fr->t0 = fr->time;
883
884 return (fr->natoms > 0);
885 }
886
887 /***** C O O R D I N A T E S T U F F *****/
888
889 int read_first_x(const output_env_t oenv,t_trxstatus **status,const char *fn,
890 real *t,rvec **x,matrix box)
891 {
892 t_trxframe fr;
893
894 read_first_frame(oenv,status,fn,&fr,TRX_NEED_X);
895
896 snew((*status)->xframe, 1);
897 (*status)->nxframe=1;
898 (*(*status)->xframe) = fr;
899 *t = (*status)->xframe->time;
900 *x = (*status)->xframe->x;
901 copy_mat((*status)->xframe->box,box);
902
903 return (*status)->xframe->natoms;
904 }
905
906 bool read_next_x(const output_env_t oenv, t_trxstatus *status,real *t,
907 int natoms, rvec x[], matrix box)
908 {
909 bool bRet;
910
911 status->xframe->x= x;
912 /*xframe[status].x = x;*/
913 bRet = read_next_frame(oenv,status,status->xframe);
914 *t = status->xframe->time;
915 copy_mat(status->xframe->box,box);
916
917 return bRet;
918 }
919
920 void close_trj(t_trxstatus *status)
921 {
922 gmx_fio_close(status->fio);
923 /* The memory in status->xframe is lost here,
924 * but the read_first_x/read_next_x functions are deprecated anyhow.
925 * read_first_frame/read_next_frame and close_trx should be used.
926 */
927 sfree(status);
928 }
929
930 void rewind_trj(t_trxstatus *status)
931 {
932 initcount(status);
933
934 gmx_fio_rewind(status->fio);
935 }
936
937 /***** V E L O C I T Y S T U F F *****/
938
939 static void clear_v(t_trxframe *fr)
940 {
941 int i;
942
943 if (!fr->bV)
944 for(i=0; i<fr->natoms; i++)
945 clear_rvec(fr->v[i]);
946 }
947
948 int read_first_v(const output_env_t oenv, t_trxstatus **status,const char *fn,
949 real *t, rvec **v,matrix box)
950 {
951 t_trxframe fr;
952
953 read_first_frame(oenv,status,fn,&fr,TRX_NEED_V);
954 *t = fr.time;
955 clear_v(&fr);
956 *v = fr.v;
957 copy_mat(fr.box,box);
958
959 return fr.natoms;
960 }
961
962 bool read_next_v(const output_env_t oenv,t_trxstatus *status,real *t,
963 int natoms,rvec v[], matrix box)
964 {
965 t_trxframe fr;
966 bool bRet;
967
968 clear_trxframe(&fr,TRUE);
969 fr.flags = TRX_NEED_V;
970 fr.natoms = natoms;
971 fr.time = *t;
972 fr.v = v;
973 bRet = read_next_frame(oenv,status,&fr);
974 *t = fr.time;
975 clear_v(&fr);
976 copy_mat(fr.box,box);
977
978 return bRet;
979 }