?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Partial commit of the project to remove all static variables.
[gromacs.git] / src / gmxlib / confio.c
blobf3c1b3ab196a6918ff9c5121e31eb77415fd6d17
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 3.1
11 * Copyright (c) 1991-2001, University of Groningen, The Netherlands
12 * This program is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU General Public License
14 * as published by the Free Software Foundation; either version 2
15 * of the License, or (at your option) any later version.
16 *
17 * If you want to redistribute modifications, please consider that
18 * scientific software is very special. Version control is crucial -
19 * bugs must be traceable. We will be happy to consider code for
20 * inclusion in the official distribution, but derived work must not
21 * be called official GROMACS. Details are found in the README & COPYING
22 * files - if they are missing, get the official version at www.gromacs.org.
23 *
24 * To help us fund GROMACS development, we humbly ask that you cite
25 * the papers on the package - you can find them in the top README file.
26 *
27 * For more info, check our website at http://www.gromacs.org
28 *
29 * And Hey:
30 * Great Red Owns Many ACres of Sand
31 */
32
33 #ifdef HAVE_CONFIG_H
34 #include <config.h>
35 #endif
36
37 #include <math.h>
38 #include "sysstuff.h"
39 #include "typedefs.h"
40 #include "smalloc.h"
41 #include "sysstuff.h"
42 #include "errno.h"
43 #include "macros.h"
44 #include "string2.h"
45 #include "confio.h"
46 #include "vec.h"
47 #include "symtab.h"
48 #include "assert.h"
49 #include "futil.h"
50 #include "xdrf.h"
51 #include "filenm.h"
52 #include "pdbio.h"
53 #include "tpxio.h"
54 #include "fatal.h"
55 #include "copyrite.h"
56 #include "filenm.h"
57 #include "statutil.h"
58
59 #define CHAR_SHIFT 24
60
61 static int read_g96_pos(char line[],t_symtab *symtab,FILE *fp,char *infile,
62 t_trxframe *fr)
63 {
64 t_atoms *atoms;
65 bool bEnd;
66 int nwanted,natoms,atnr,resnr,oldres,newres,shift;
67 char anm[STRLEN],resnm[STRLEN];
68 char c1,c2;
69 double db1,db2,db3;
70
71 nwanted = fr->natoms;
72
73 atoms = fr->atoms;
74
75 natoms = 0;
76
77 if (fr->bX) {
78 if (fr->bAtoms)
79 shift = CHAR_SHIFT;
80 else
81 shift = 0;
82 newres = 0;
83 oldres = -666; /* Unlikely number for the first residue! */
84 bEnd = FALSE;
85 while (!bEnd && fgets2(line,STRLEN,fp)) {
86 bEnd = (strncmp(line,"END",3) == 0);
87 if (!bEnd && (line[0] != '#')) {
88 if (sscanf(line+shift,"%15lf%15lf%15lf",&db1,&db2,&db3) != 3)
89 fatal_error(0,"Did not find 3 coordinates for atom %d in %s\n",
90 natoms+1,infile);
91 if ((nwanted != -1) && (natoms >= nwanted))
92 fatal_error(0,
93 "Found more coordinates (%d) in %s than expected %d\n",
94 natoms,infile,nwanted);
95 if (atoms) {
96 if (atoms && fr->bAtoms &&
97 (sscanf(line,"%5d%c%5s%c%5s%7d",&resnr,&c1,resnm,&c2,anm,&atnr)
98 != 6)) {
99 if (oldres>=0)
100 resnr = oldres;
101 else {
102 resnr = 1;
103 strcpy(resnm,"???");
104 }
105 strcpy(anm,"???");
106 }
107 atoms->atomname[natoms]=put_symtab(symtab,anm);
108 if (resnr != oldres) {
109 oldres = resnr;
110 if (newres >= atoms->nr)
111 fatal_error(0,"More residues than atoms in %s (natoms = %d)",
112 infile,atoms->nr);
113 atoms->resname[newres] = put_symtab(symtab,resnm);
114 newres++;
115 if (newres > atoms->nres)
116 atoms->nres = newres;
117 }
118 resnr = newres;
119 atoms->atom[natoms].resnr = resnr-1;
120 }
121 if (fr->x) {
122 fr->x[natoms][0] = db1;
123 fr->x[natoms][1] = db2;
124 fr->x[natoms][2] = db3;
125 }
126 natoms++;
127 }
128 }
129 if ((nwanted != -1) && natoms != nwanted)
130 fprintf(stderr,
131 "Warning: found less coordinates (%d) in %s than expected %d\n",
132 natoms,infile,nwanted);
133 }
134
135 fr->natoms = natoms;
136
137 return natoms;
138 }
139
140 static int read_g96_vel(char line[],FILE *fp,char *infile,
141 t_trxframe *fr)
142 {
143 bool bEnd;
144 int nwanted,natoms=-1,shift;
145 double db1,db2,db3;
146
147 nwanted = fr->natoms;
148
149 if (fr->v && fr->bV) {
150 if (strcmp(line,"VELOCITYRED") == 0)
151 shift = 0;
152 else
153 shift = CHAR_SHIFT;
154 natoms = 0;
155 bEnd = FALSE;
156 while (!bEnd && fgets2(line,STRLEN,fp)) {
157 bEnd = (strncmp(line,"END",3) == 0);
158 if (!bEnd && (line[0] != '#')) {
159 if (sscanf(line+shift,"%15lf%15lf%15lf",&db1,&db2,&db3) != 3)
160 fatal_error(0,"Did not find 3 velocities for atom %d in %s",
161 natoms+1,infile);
162 if ((nwanted != -1) && (natoms >= nwanted))
163 fatal_error(0,"Found more velocities (%d) in %s than expected %d\n",
164 natoms,infile,nwanted);
165 if (fr->v) {
166 fr->v[natoms][0] = db1;
167 fr->v[natoms][1] = db2;
168 fr->v[natoms][2] = db3;
169 }
170 natoms++;
171 }
172 }
173 if ((nwanted != -1) && (natoms != nwanted))
174 fprintf(stderr,
175 "Warning: found less velocities (%d) in %s than expected %d\n",
176 natoms,infile,nwanted);
177 }
178
179 return natoms;
180 }
181
182 int read_g96_conf(FILE *fp,char *infile,t_trxframe *fr)
183 {
184 static t_symtab *symtab=NULL;
185 static char line[STRLEN+1]; /* VERY DIRTY, you can not read two *
186 * Gromos96 trajectories at the same time */
187 bool bAtStart,bTime,bAtoms,bPos,bVel,bBox,bEnd,bFinished;
188 int natoms,nbp;
189 double db1,db2,db3,db4,db5,db6,db7,db8,db9;
190
191 bAtStart = (ftell(fp) == 0);
192
193 clear_trxframe(fr,FALSE);
194
195 if (!symtab) {
196 snew(symtab,1);
197 open_symtab(symtab);
198 }
199
200 natoms=0;
201
202 if (bAtStart) {
203 while ( !fr->bTitle && fgets2(line,STRLEN,fp))
204 fr->bTitle = (strcmp(line,"TITLE") == 0);
205 if (fr->title)
206 fgets2(fr->title,STRLEN,fp);
207 bEnd = FALSE;
208 while (!bEnd && fgets2(line,STRLEN,fp))
209 bEnd = (strcmp(line,"END") == 0);
210 fgets2(line,STRLEN,fp);
211 }
212
213 /* Do not get a line if we are not at the start of the file, *
214 * because without a parameter file we don't know what is in *
215 * the trajectory and we have already read the line in the *
216 * previous call (VERY DIRTY). */
217 bFinished = FALSE;
218 do {
219 bTime = (strcmp(line,"TIMESTEP") == 0);
220 bAtoms = (strcmp(line,"POSITION") == 0);
221 bPos = (bAtoms || (strcmp(line,"POSITIONRED") == 0));
222 bVel = (strncmp(line,"VELOCITY",8) == 0);
223 bBox = (strcmp(line,"BOX") == 0);
224 if (bTime) {
225 if (!fr->bTime && !fr->bX) {
226 fr->bStep = bTime;
227 fr->bTime = bTime;
228 do
229 bFinished = (fgets2(line,STRLEN,fp) == NULL);
230 while (!bFinished && (line[0] == '#'));
231 sscanf(line,"%15d%15lf",&(fr->step),&db1);
232 fr->time = db1;
233 } else
234 bFinished = TRUE;
235 }
236 if (bPos) {
237 if (!fr->bX) {
238 fr->bAtoms = bAtoms;
239 fr->bX = bPos;
240 natoms = read_g96_pos(line,symtab,fp,infile,fr);
241 } else
242 bFinished = TRUE;
243 }
244 if (fr->v && bVel) {
245 fr->bV = bVel;
246 natoms = read_g96_vel(line,fp,infile,fr);
247 }
248 if (bBox) {
249 fr->bBox = bBox;
250 clear_mat(fr->box);
251 bEnd = FALSE;
252 while (!bEnd && fgets2(line,STRLEN,fp)) {
253 bEnd = (strncmp(line,"END",3) == 0);
254 if (!bEnd && (line[0] != '#')) {
255 nbp = sscanf(line,"%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf",
256 &db1,&db2,&db3,&db4,&db5,&db6,&db7,&db8,&db9);
257 if (nbp < 3)
258 fatal_error(0,"Found a BOX line, but no box in %s",infile);
259 fr->box[XX][XX] = db1;
260 fr->box[YY][YY] = db2;
261 fr->box[ZZ][ZZ] = db3;
262 if (nbp == 9) {
263 fr->box[XX][YY] = db4;
264 fr->box[XX][ZZ] = db5;
265 fr->box[YY][XX] = db6;
266 fr->box[YY][ZZ] = db7;
267 fr->box[ZZ][XX] = db8;
268 fr->box[ZZ][YY] = db9;
269 }
270 }
271 }
272 bFinished = TRUE;
273 }
274 } while (!bFinished && fgets2(line,STRLEN,fp));
275
276 close_symtab(symtab);
277
278 fr->natoms = natoms;
279
280 return natoms;
281 }
282
283 void write_g96_conf(FILE *out,t_trxframe *fr,
284 int nindex,atom_id *index)
285 {
286 t_atoms *atoms;
287 int nout,i,a;
288
289 atoms = fr->atoms;
290
291 if (index)
292 nout = nindex;
293 else
294 nout = fr->natoms;
295
296 if (fr->bTitle)
297 fprintf(out,"TITLE\n%s\nEND\n",fr->title);
298 if (fr->bStep || fr->bTime)
299 /* Officially the time format is %15.9, which is not enough for 10 ns */
300 fprintf(out,"TIMESTEP\n%15d%15.6f\nEND\n",fr->step,fr->time);
301 if (fr->bX) {
302 if (fr->bAtoms) {
303 fprintf(out,"POSITION\n");
304 for(i=0; i<nout; i++) {
305 if (index) a = index[i]; else a = i;
306 fprintf(out,"%5d %-5s %-5s%7d%15.9f%15.9f%15.9f\n",
307 (atoms->atom[a].resnr+1) % 100000,
308 *atoms->resname[atoms->atom[a].resnr],
309 *atoms->atomname[a],(i+1) % 10000000,
310 fr->x[a][XX],fr->x[a][YY],fr->x[a][ZZ]);
311 }
312 } else {
313 fprintf(out,"POSITIONRED\n");
314 for(i=0; i<nout; i++) {
315 if (index) a = index[i]; else a = i;
316 fprintf(out,"%15.9f%15.9f%15.9f\n",
317 fr->x[a][XX],fr->x[a][YY],fr->x[a][ZZ]);
318 }
319 }
320 fprintf(out,"END\n");
321 }
322 if (fr->bV) {
323 if (atoms) {
324 fprintf(out,"VELOCITY\n");
325 for(i=0; i<nout; i++) {
326 if (index) a = index[i]; else a = i;
327 fprintf(out,"%5d %-5s %-5s%7d%15.9f%15.9f%15.9f\n",
328 (atoms->atom[a].resnr+1) % 100000,
329 *atoms->resname[atoms->atom[a].resnr],
330 *atoms->atomname[a],(i+1) % 10000000,
331 fr->v[a][XX],fr->v[a][YY],fr->v[a][ZZ]);
332 }
333 } else {
334 fprintf(out,"VELOCITYRED\n");
335 for(i=0; i<nout; i++) {
336 if (index) a = index[i]; else a = i;
337 fprintf(out,"%15.9f%15.9f%15.9f\n",
338 fr->v[a][XX],fr->v[a][YY],fr->v[a][ZZ]);
339 }
340 }
341 fprintf(out,"END\n");
342 }
343 if (fr->bBox) {
344 fprintf(out,"BOX\n");
345 fprintf(out,"%15.9f%15.9f%15.9f",
346 fr->box[XX][XX],fr->box[YY][YY],fr->box[ZZ][ZZ]);
347 if (fr->box[XX][YY] || fr->box[XX][ZZ] || fr->box[YY][XX] ||
348 fr->box[YY][ZZ] || fr->box[ZZ][XX] ||fr->box[ZZ][YY])
349 fprintf(out,"%15.9f%15.9f%15.9f%15.9f%15.9f%15.9f",
350 fr->box[XX][YY],fr->box[XX][ZZ],fr->box[YY][XX],
351 fr->box[YY][ZZ],fr->box[ZZ][XX],fr->box[ZZ][YY]);
352 fprintf(out,"\n");
353 fprintf(out,"END\n");
354 }
355 }
356
357 static void get_coordnum_fp (FILE *in,char *title, int *natoms)
358 {
359 char line[STRLEN+1];
360
361 fgets2 (title,STRLEN,in);
362 fgets2 (line,STRLEN,in);
363 sscanf (line,"%d",natoms);
364 }
365
366 static void get_coordnum (char *infile,int *natoms)
367 {
368 FILE *in;
369 char title[STRLEN];
370
371 in=ffopen(infile,"r");
372 get_coordnum_fp(in,title,natoms);
373 ffclose (in);
374 }
375
376 static bool get_w_conf(FILE *in, char *infile, char *title,
377 t_atoms *atoms, int *ndec, rvec x[],rvec *v, matrix box)
378 {
379 static t_symtab *symtab=NULL;
380 char name[6];
381 char line[STRLEN+1];
382 char buf[256];
383 char format[30];
384 double x1,y1,z1,x2,y2,z2;
385 rvec xmin,xmax;
386 int natoms,i,m,resnr,newres,oldres,ddist;
387 bool bFirst,bVel;
388 char *p1,*p2;
389
390 newres = 0;
391 oldres = NOTSET; /* Unlikely number for the first residue! */
392 ddist = 0;
393
394 if (!symtab) {
395 snew(symtab,1);
396 open_symtab(symtab);
397 }
398 fgets2(title,STRLEN,in);
399
400 /* read the number of atoms */
401 fgets2(line,STRLEN,in);
402 sscanf(line,"%d",&natoms);
403 if (natoms > atoms->nr)
404 fatal_error(0,"gro file contains more atoms (%d) than expected (%d)",
405 natoms,atoms->nr);
406 else if (natoms < atoms->nr)
407 fprintf(stderr,"Warning: gro file contains less atoms (%d) than expected"
408 " (%d)\n",natoms,atoms->nr);
409
410 bFirst=TRUE;
411
412 bVel = FALSE;
413
414 /* just pray the arrays are big enough */
415 for (i=0; (i < natoms) ; i++) {
416 if ((fgets2 (line,STRLEN,in)) == NULL) {
417 unexpected_eof(infile,i+2);
418 }
419 if (strlen(line) < 39)
420 fatal_error(0,"Invalid line in %s for atom %d:\n%s",infile,i+1,line);
421
422 /* determine read precision from distance between periods
423 (decimal points) */
424 if (bFirst) {
425 bFirst=FALSE;
426 p1=strchr(line,'.');
427 if (p1 == NULL)
428 fatal_error(0,"A coordinate in file %s does not contain a '.'",infile);
429 p2=strchr(&p1[1],'.');
430 if (p1 || p2)
431 ddist=p2-p1;
432 else
433 ddist=8;
434 if (ddist<0)
435 ddist=8;
436 if (ddist>30)
437 ddist=30;
438 sprintf(format,"%%%dlf%%%dlf%%%dlf",ddist,ddist,ddist);
439 /* this will be something like "%8lf%8lf%8lf" */
440 *ndec = ddist-5;
441 }
442
443 /* residue number*/
444 memcpy(name,line,5);
445 name[5]='\0';
446 sscanf(name,"%d",&resnr);
447 memcpy(name,line+5,5);
448 name[5]='\0';
449 if (resnr != oldres) {
450 oldres = resnr;
451 if (newres >= natoms)
452 fatal_error(0,"More residues than atoms in %s (natoms = %d)",
453 infile,natoms);
454 atoms->resname[newres] = put_symtab(symtab,name);
455 newres++;
456 }
457 resnr = newres;
458 atoms->atom[i].resnr = resnr-1;
459
460 /* atomname */
461 memcpy(name,line+10,5);
462 atoms->atomname[i]=put_symtab(symtab,name);
463
464 /* eventueel controle atomnumber met i+1 */
465
466 /* coordinates (start after residue shit) */
467 /* 'format' was built previously */
468 if (sscanf (line+20,format,&x1,&y1,&z1) != 3) {
469 too_few();
470 }
471 else {
472 x[i][XX]=x1;
473 x[i][YY]=y1;
474 x[i][ZZ]=z1;
475 }
476
477 /* velocities (start after residues and coordinates) */
478 /* 'format' was built previously */
479 if (v) {
480 if (sscanf (line+20+(3*ddist),format,&x1,&y1,&z1) != 3) {
481 v[i][XX] = 0.0;
482 v[i][YY] = 0.0;
483 v[i][ZZ] = 0.0;
484 }
485 else {
486 v[i][XX]=x1;
487 v[i][YY]=y1;
488 v[i][ZZ]=z1;
489 bVel = TRUE;
490 }
491 }
492 }
493 atoms->nres=newres;
494
495 /* box */
496 fgets2 (line,STRLEN,in);
497 if (sscanf (line,"%lf%lf%lf",&x1,&y1,&z1) != 3) {
498 sprintf(buf,"Bad box in file %s",infile);
499 warning(buf);
500
501 /* Generate a cubic box */
502 for(m=0; (m<DIM); m++)
503 xmin[m]=xmax[m]=x[0][m];
504 for(i=1; (i<atoms->nr); i++)
505 for(m=0; (m<DIM); m++) {
506 xmin[m]=min(xmin[m],x[i][m]);
507 xmax[m]=max(xmax[m],x[i][m]);
508 }
509 for (i=0; i<DIM; i++) for (m=0; m<DIM; m++) box[i][m]=0.0;
510 for(m=0; (m<DIM); m++)
511 box[m][m]=(xmax[m]-xmin[m]);
512 fprintf(stderr,"Generated a cubic box %8.3f x %8.3f x %8.3f\n",
513 box[XX][XX],box[YY][YY],box[ZZ][ZZ]);
514 }
515 else {
516 /* We found the first three values, the diagonal elements */
517 box[XX][XX]=x1;
518 box[YY][YY]=y1;
519 box[ZZ][ZZ]=z1;
520 if (sscanf (line,"%*f%*f%*f%lf%lf%lf%lf%lf%lf",
521 &x1,&y1,&z1,&x2,&y2,&z2) != 6)
522 x1=y1=z1=x2=y2=z2=0.0;
523 box[XX][YY] = x1;
524 box[XX][ZZ] = y1;
525 box[YY][XX] = z1;
526 box[YY][ZZ] = x2;
527 box[ZZ][XX] = y2;
528 box[ZZ][YY] = z2;
529 }
530 close_symtab(symtab);
531
532 return bVel;
533 }
534
535 static void read_whole_conf(char *infile, char *title,
536 t_atoms *atoms, rvec x[],rvec *v, matrix box)
537 {
538 FILE *in;
539 int ndec;
540
541 /* open file */
542 in=ffopen(infile,"r");
543
544 get_w_conf(in, infile, title, atoms, &ndec, x, v, box);
545
546 fclose(in);
547 }
548
549 static void get_conf(FILE *in, char *title, int *natoms,
550 rvec x[],rvec *v,matrix box)
551 {
552 t_atoms atoms;
553 int ndec;
554
555 atoms.nr=*natoms;
556 snew(atoms.atom,*natoms);
557 atoms.nres=*natoms;
558 snew(atoms.resname,*natoms);
559 snew(atoms.atomname,*natoms);
560
561 get_w_conf(in,title,title,&atoms,&ndec,x,v,box);
562
563 sfree(atoms.atom);
564 sfree(atoms.resname);
565 sfree(atoms.atomname);
566 }
567
568 bool gro_next_x_or_v(FILE *status,t_trxframe *fr)
569 {
570 t_atoms atoms;
571 char title[STRLEN],*p;
572 double tt;
573 int ndec,i;
574
575 if (eof(status))
576 return FALSE;
577
578 atoms.nr=fr->natoms;
579 snew(atoms.atom,fr->natoms);
580 atoms.nres=fr->natoms;
581 snew(atoms.resname,fr->natoms);
582 snew(atoms.atomname,fr->natoms);
583
584 fr->bV = get_w_conf(status,title,title,&atoms,&ndec,fr->x,fr->v,fr->box);
585 fr->bPrec = TRUE;
586 fr->prec = 1;
587 /* prec = 10^ndec: */
588 for(i=0; i<ndec; i++)
589 fr->prec *= 10;
590 fr->title = title;
591 fr->bTitle = TRUE;
592 fr->bX = TRUE;
593 fr->bBox = TRUE;
594
595 sfree(atoms.atom);
596 sfree(atoms.resname);
597 sfree(atoms.atomname);
598
599 if ((p=strstr(title,"t=")) != NULL) {
600 p+=2;
601 if (sscanf(p,"%lf",&tt)==1) {
602 fr->time = tt;
603 fr->bTime = TRUE;
604 } else {
605 fr->time = 0;
606 fr->bTime = FALSE;
607 }
608 }
609
610 if (atoms.nr != fr->natoms)
611 fatal_error(0,"Number of atoms in gro frame (%d) doesn't match the number in the previous frame (%d)",atoms.nr,fr->natoms);
612
613 return TRUE;
614 }
615
616 int gro_first_x_or_v(FILE *status,t_trxframe *fr)
617 {
618 char title[STRLEN];
619
620 frewind(status);
621 fprintf(stderr,"Reading frames from gro file");
622 get_coordnum_fp(status, title, &fr->natoms);
623 frewind(status);
624 fprintf(stderr," '%s', %d atoms.\n",title, fr->natoms);
625 fr->bTitle = TRUE;
626 fr->title = title;
627 if (fr->natoms==0)
628 fatal_error(1,"No coordinates in gro file\n");
629
630 snew(fr->x,fr->natoms);
631 snew(fr->v,fr->natoms);
632 gro_next_x_or_v(status, fr);
633
634 return fr->natoms;
635 }
636
637 void write_hconf_indexed_p(FILE *out,char *title,t_atoms *atoms,
638 int nx,atom_id index[], int pr,
639 rvec *x,rvec *v,matrix box)
640 {
641 char resnm[6],nm[6],format[100];
642 int ai,i,resnr,l,vpr;
643
644 fprintf (out,"%s\n",(title && title[0])?title:bromacs(format,99));
645 fprintf (out,"%5d\n",nx);
646 /* build format string for printing,
647 something like "%8.3f" for x and "%8.4f" for v */
648 if (pr<0)
649 pr=0;
650 if (pr>30)
651 pr=30;
652 l=pr+5;
653 vpr=pr+1;
654 if (v)
655 sprintf(format,"%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df\n",
656 l,pr,l,pr,l,pr,l,vpr,l,vpr,l,vpr);
657 else
658 sprintf(format,"%%%d.%df%%%d.%df%%%d.%df\n",l,pr,l,pr,l,pr);
659
660 for (i=0; (i<nx); i++) {
661 ai=index[i];
662
663 resnr=atoms->atom[ai].resnr;
664 strcpy(resnm," ??? ");
665 if (resnr < atoms->nres)
666 strcpy(resnm,*atoms->resname[resnr]);
667
668 if (atoms->atom)
669 strcpy(nm,*atoms->atomname[ai]);
670 else
671 strcpy(nm," ??? ");
672
673 fprintf(out,"%5d%-5.5s%5.5s%5d",(resnr+1)%100000,resnm,nm,(ai+1)%100000);
674 /* next fprintf uses built format string */
675 if (v)
676 fprintf(out,format,
677 x[ai][XX], x[ai][YY], x[ai][ZZ], v[ai][XX],v[ai][YY],v[ai][ZZ]);
678 else
679 fprintf(out,format,
680 x[ai][XX], x[ai][YY], x[ai][ZZ]);
681 }
682
683 if (pr<5)
684 pr=5;
685 l=pr+5;
686
687 if (box[XX][YY] || box[XX][ZZ] || box[YY][XX] || box[YY][ZZ] ||
688 box[ZZ][XX] || box[ZZ][YY]) {
689 sprintf(format,"%%%d.%df%%%d.%df%%%d.%df"
690 "%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df\n",
691 l,pr,l,pr,l,pr,l,pr,l,pr,l,pr,l,pr,l,pr,l,pr);
692 fprintf(out,format,
693 box[XX][XX],box[YY][YY],box[ZZ][ZZ],
694 box[XX][YY],box[XX][ZZ],box[YY][XX],
695 box[YY][ZZ],box[ZZ][XX],box[ZZ][YY]);
696 } else {
697 sprintf(format,"%%%d.%df%%%d.%df%%%d.%df\n",l,pr,l,pr,l,pr);
698 fprintf(out,format,
699 box[XX][XX],box[YY][YY],box[ZZ][ZZ]);
700 }
701 fflush(out);
702 }
703
704 void write_hconf_p(FILE *out,char *title,t_atoms *atoms, int pr,
705 rvec *x,rvec *v,matrix box)
706 {
707 atom_id *aa;
708 int i;
709
710 snew(aa,atoms->nr);
711 for(i=0; (i<atoms->nr); i++)
712 aa[i]=i;
713 write_hconf_indexed_p(out,title,atoms,atoms->nr,aa,pr,x,v,box);
714 sfree(aa);
715 }
716
717 void write_conf_p(char *outfile, char *title, t_atoms *atoms, int pr,
718 rvec *x, rvec *v,matrix box)
719 {
720 FILE *out;
721
722 out=ffopen(outfile,"w");
723 write_hconf_p(out,title,atoms,pr,x,v,box);
724
725 ffclose (out);
726 }
727
728 static void write_conf(char *outfile, char *title, t_atoms *atoms,
729 rvec *x, rvec *v,matrix box)
730 {
731 write_conf_p(outfile, title, atoms, 3, x, v, box);
732 }
733
734 void write_sto_conf_indexed(char *outfile,char *title,t_atoms *atoms,
735 rvec x[],rvec *v,matrix box,
736 atom_id nindex,atom_id index[])
737 {
738 FILE *out;
739 int ftp;
740 t_trxframe fr;
741
742 ftp=fn2ftp(outfile);
743 switch (ftp) {
744 case efGRO:
745 out=ffopen(outfile,"w");
746 write_hconf_indexed_p(out, title, atoms, nindex, index, 3, x, v, box);
747 fclose(out);
748 break;
749 case efG96:
750 clear_trxframe(&fr,TRUE);
751 fr.bTitle = TRUE;
752 fr.title = title;
753 fr.natoms = atoms->nr;
754 fr.bAtoms = TRUE;
755 fr.atoms = atoms;
756 fr.bX = TRUE;
757 fr.x = x;
758 if (v) {
759 fr.bV = TRUE;
760 fr.v = v;
761 }
762 fr.bBox = TRUE;
763 copy_mat(box,fr.box);
764 out=ffopen(outfile,"w");
765 write_g96_conf(out, &fr, nindex, index);
766 fclose(out);
767 break;
768 case efPDB:
769 case efBRK:
770 case efENT:
771 out=ffopen(outfile,"w");
772 write_pdbfile_indexed(out, title, atoms, x, box, 0, -1, nindex, index);
773 fclose(out);
774 break;
775 case efTPR:
776 case efTPB:
777 case efTPA:
778 fatal_error(0,"Sorry, can not write a topology to %s",outfile);
779 break;
780 default:
781 fatal_error(0,"Not supported in write_sto_conf_indexed: %s",outfile);
782 }
783 }
784
785 void write_sto_conf(char *outfile, char *title,t_atoms *atoms,
786 rvec x[],rvec *v, matrix box)
787 {
788 FILE *out;
789 int ftp;
790 t_trxframe fr;
791
792 ftp=fn2ftp(outfile);
793 switch (ftp) {
794 case efGRO:
795 write_conf(outfile, title, atoms, x, v, box);
796 break;
797 case efG96:
798 clear_trxframe(&fr,TRUE);
799 fr.bTitle = TRUE;
800 fr.title = title;
801 fr.natoms = atoms->nr;
802 fr.bAtoms = TRUE;
803 fr.atoms = atoms;
804 fr.bX = TRUE;
805 fr.x = x;
806 if (v) {
807 fr.bV = TRUE;
808 fr.v = v;
809 }
810 fr.bBox = TRUE;
811 copy_mat(box,fr.box);
812 out=ffopen(outfile,"w");
813 write_g96_conf(out, &fr, -1, NULL);
814 fclose(out);
815 break;
816 case efPDB:
817 case efBRK:
818 case efENT:
819 out=ffopen(outfile,"w");
820 write_pdbfile(out, title, atoms, x, box, 0, -1);
821 fclose(out);
822 break;
823 case efTPR:
824 case efTPB:
825 case efTPA:
826 fatal_error(0,"Sorry, can not write a topology to %s",outfile);
827 break;
828 default:
829 fatal_error(0,"Not supported in write_sto_conf: %s",outfile);
830 }
831 }
832
833 void get_stx_coordnum(char *infile,int *natoms)
834 {
835 FILE *in;
836 int ftp,tpxver,tpxgen;
837 t_trxframe fr;
838
839 ftp=fn2ftp(infile);
840 switch (ftp) {
841 case efGRO:
842 get_coordnum(infile, natoms);
843 break;
844 case efG96:
845 in=ffopen(infile,"r");
846 fr.title = NULL;
847 fr.natoms = -1;
848 fr.atoms = NULL;
849 fr.x = NULL;
850 fr.v = NULL;
851 fr.f = NULL;
852 *natoms=read_g96_conf(in,infile,&fr);
853 fclose(in);
854 break;
855 case efPDB:
856 case efBRK:
857 case efENT:
858 in=ffopen(infile,"r");
859 get_pdb_coordnum(in, natoms);
860 fclose(in);
861 break;
862 case efTPA:
863 case efTPB:
864 case efTPR: {
865 t_tpxheader tpx;
866
867 read_tpxheader(infile,&tpx,TRUE,&tpxver,&tpxgen);
868 *natoms = tpx.natoms;
869 break;
870 }
871 default:
872 fatal_error(0,"Not supported in get_stx_coordnum: %s",infile);
873 }
874 }
875
876 void read_stx_conf(char *infile, char *title,t_atoms *atoms,
877 rvec x[],rvec *v, matrix box)
878 {
879 FILE *in;
880 char buf[256];
881 t_topology *top;
882 t_trxframe fr;
883 int i,ftp,natoms,i1;
884 real d,r1,r2;
885
886 if (atoms->nr == 0)
887 fprintf(stderr,"Warning: Number of atoms in %s is 0\n",infile);
888 else if (atoms->atom == NULL) {
889 sprintf(buf,"Uninitialized array atom in %s, %d",__FILE__,__LINE__);
890 fatal_error(0,buf);
891 }
892 ftp=fn2ftp(infile);
893 switch (ftp) {
894 case efGRO:
895 read_whole_conf(infile, title, atoms, x, v, box);
896 break;
897 case efG96:
898 fr.title = title;
899 fr.natoms = atoms->nr;
900 fr.atoms = atoms;
901 fr.x = x;
902 fr.v = v;
903 fr.f = NULL;
904 in = ffopen(infile,"r");
905 read_g96_conf(in, infile, &fr);
906 fclose(in);
907 copy_mat(fr.box,box);
908 break;
909 case efPDB:
910 case efBRK:
911 case efENT:
912 read_pdb_conf(infile, title, atoms, x, box, TRUE);
913 break;
914 case efTPR:
915 case efTPB:
916 case efTPA:
917 snew(top,1);
918 read_tpx(infile,&i1,&r1,&r2,NULL,box,&natoms,x,v,NULL,top);
919
920 strcpy(title,*(top->name));
921 /* Scalars */
922 atoms->nr = top->atoms.nr;
923 atoms->nres = top->atoms.nres;
924 atoms->ngrpname = top->atoms.ngrpname;
925
926 /* Arrays */
927 if (!atoms->atom)
928 snew(atoms->atom,atoms->nr);
929 if (!atoms->atomname)
930 snew(atoms->atomname,atoms->nr);
931 if (!atoms->atomtype)
932 snew(atoms->atomtype,atoms->nr);
933 if (!atoms->atomtypeB)
934 snew(atoms->atomtypeB,atoms->nr);
935 for(i=0; (i<atoms->nr); i++) {
936 atoms->atom[i] = top->atoms.atom[i];
937 atoms->atomname[i] = top->atoms.atomname[i];
938 atoms->atomtype[i] = top->atoms.atomtype[i];
939 atoms->atomtypeB[i] = top->atoms.atomtypeB[i];
940
941 }
942
943 if (!atoms->resname)
944 snew(atoms->resname,atoms->nres);
945 for(i=0; (i<atoms->nres); i++)
946 atoms->resname[i] = top->atoms.resname[i];
947
948 if (!atoms->grpname)
949 snew(atoms->grpname,atoms->ngrpname);
950 for(i=0; (i<atoms->ngrpname); i++)
951 atoms->grpname[i] = top->atoms.grpname[i];
952
953 for(i=0; (i<egcNR); i++) {
954 atoms->grps[i].nr = top->atoms.grps[i].nr;
955 if (atoms->grps[i].nr > 0) {
956 snew(atoms->grps[i].nm_ind,atoms->grps[i].nr);
957 memcpy(atoms->grps[i].nm_ind,top->atoms.grps[i].nm_ind,
958 atoms->grps[i].nr*sizeof(atoms->grps[i].nm_ind[0]));
959 }
960 }
961
962 /* Ignore exclusions */
963
964 break;
965 default:
966 fatal_error(0,"Not supported in read_stx_conf: %s",infile);
967 }
968 }
969
970
The ultimate molecular dynamics simulation package