?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Revised wording in pdb2gmx.c, hopefully clearer now.
[gromacs/rigid-bodies.git] / src / gmxlib / confio.c
blob2acf638dbd89d07da7dbb0e7095d8ea912086f36
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 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.
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 <math.h>
40 #include "sysstuff.h"
41 #include "typedefs.h"
42 #include "smalloc.h"
43 #include "sysstuff.h"
44 #include "errno.h"
45 #include "macros.h"
46 #include "string2.h"
47 #include "confio.h"
48 #include "vec.h"
49 #include "symtab.h"
50 #include "futil.h"
51 #include "xdrf.h"
52 #include "filenm.h"
53 #include "pdbio.h"
54 #include "tpxio.h"
55 #include "gmx_fatal.h"
56 #include "copyrite.h"
57 #include "filenm.h"
58 #include "statutil.h"
59 #include "pbc.h"
60 #include "mtop_util.h"
61 #include "gmxfio.h"
62
63 #define CHAR_SHIFT 24
64
65 static int read_g96_pos(char line[],t_symtab *symtab,
66 FILE *fp,const char *infile,
67 t_trxframe *fr)
68 {
69 t_atoms *atoms;
70 gmx_bool bEnd;
71 int nwanted,natoms,atnr,resnr,oldres,newres,shift;
72 char anm[STRLEN],resnm[STRLEN];
73 char c1,c2;
74 double db1,db2,db3;
75
76 nwanted = fr->natoms;
77
78 atoms = fr->atoms;
79
80 natoms = 0;
81
82 if (fr->bX) {
83 if (fr->bAtoms)
84 shift = CHAR_SHIFT;
85 else
86 shift = 0;
87 newres = -1;
88 oldres = -666; /* Unlikely number for the first residue! */
89 bEnd = FALSE;
90 while (!bEnd && fgets2(line,STRLEN,fp)) {
91 bEnd = (strncmp(line,"END",3) == 0);
92 if (!bEnd && (line[0] != '#')) {
93 if (sscanf(line+shift,"%15lf%15lf%15lf",&db1,&db2,&db3) != 3)
94 gmx_fatal(FARGS,"Did not find 3 coordinates for atom %d in %s\n",
95 natoms+1,infile);
96 if ((nwanted != -1) && (natoms >= nwanted))
97 gmx_fatal(FARGS,
98 "Found more coordinates (%d) in %s than expected %d\n",
99 natoms,infile,nwanted);
100 if (atoms) {
101 if (atoms && fr->bAtoms &&
102 (sscanf(line,"%5d%c%5s%c%5s%7d",&resnr,&c1,resnm,&c2,anm,&atnr)
103 != 6)) {
104 if (oldres>=0)
105 resnr = oldres;
106 else {
107 resnr = 1;
108 strncpy(resnm,"???",sizeof(resnm)-1);
109 }
110 strncpy(anm,"???",sizeof(anm)-1);
111 }
112 atoms->atomname[natoms]=put_symtab(symtab,anm);
113 if (resnr != oldres) {
114 oldres = resnr;
115 newres++;
116 if (newres >= atoms->nr)
117 gmx_fatal(FARGS,"More residues than atoms in %s (natoms = %d)",
118 infile,atoms->nr);
119 atoms->atom[natoms].resind = newres;
120 if (newres+1 > atoms->nres) {
121 atoms->nres = newres+1;
122 }
123 t_atoms_set_resinfo(atoms,natoms,symtab,resnm,resnr,' ',0,' ');
124 } else {
125 atoms->atom[natoms].resind = newres;
126 }
127 }
128 if (fr->x) {
129 fr->x[natoms][0] = db1;
130 fr->x[natoms][1] = db2;
131 fr->x[natoms][2] = db3;
132 }
133 natoms++;
134 }
135 }
136 if ((nwanted != -1) && natoms != nwanted)
137 fprintf(stderr,
138 "Warning: found less coordinates (%d) in %s than expected %d\n",
139 natoms,infile,nwanted);
140 }
141
142 fr->natoms = natoms;
143
144 return natoms;
145 }
146
147 static int read_g96_vel(char line[],FILE *fp,const char *infile,
148 t_trxframe *fr)
149 {
150 gmx_bool bEnd;
151 int nwanted,natoms=-1,shift;
152 double db1,db2,db3;
153
154 nwanted = fr->natoms;
155
156 if (fr->v && fr->bV) {
157 if (strcmp(line,"VELOCITYRED") == 0)
158 shift = 0;
159 else
160 shift = CHAR_SHIFT;
161 natoms = 0;
162 bEnd = FALSE;
163 while (!bEnd && fgets2(line,STRLEN,fp)) {
164 bEnd = (strncmp(line,"END",3) == 0);
165 if (!bEnd && (line[0] != '#')) {
166 if (sscanf(line+shift,"%15lf%15lf%15lf",&db1,&db2,&db3) != 3)
167 gmx_fatal(FARGS,"Did not find 3 velocities for atom %d in %s",
168 natoms+1,infile);
169 if ((nwanted != -1) && (natoms >= nwanted))
170 gmx_fatal(FARGS,"Found more velocities (%d) in %s than expected %d\n",
171 natoms,infile,nwanted);
172 if (fr->v) {
173 fr->v[natoms][0] = db1;
174 fr->v[natoms][1] = db2;
175 fr->v[natoms][2] = db3;
176 }
177 natoms++;
178 }
179 }
180 if ((nwanted != -1) && (natoms != nwanted))
181 fprintf(stderr,
182 "Warning: found less velocities (%d) in %s than expected %d\n",
183 natoms,infile,nwanted);
184 }
185
186 return natoms;
187 }
188
189 int read_g96_conf(FILE *fp,const char *infile,t_trxframe *fr)
190 {
191 t_symtab *symtab=NULL;
192 char line[STRLEN+1];
193 gmx_bool bAtStart,bTime,bAtoms,bPos,bVel,bBox,bEnd,bFinished;
194 int natoms,nbp;
195 double db1,db2,db3,db4,db5,db6,db7,db8,db9;
196
197 bAtStart = (ftell(fp) == 0);
198
199 clear_trxframe(fr,FALSE);
200
201 if (!symtab) {
202 snew(symtab,1);
203 open_symtab(symtab);
204 }
205
206 natoms=0;
207
208 if (bAtStart) {
209 while ( !fr->bTitle && fgets2(line,STRLEN,fp))
210 fr->bTitle = (strcmp(line,"TITLE") == 0);
211 if (fr->title == NULL) {
212 fgets2(line,STRLEN,fp);
213 fr->title = strdup(line);
214 }
215 bEnd = FALSE;
216 while (!bEnd && fgets2(line,STRLEN,fp))
217 bEnd = (strcmp(line,"END") == 0);
218 fgets2(line,STRLEN,fp);
219 }
220
221 /* Do not get a line if we are not at the start of the file, *
222 * because without a parameter file we don't know what is in *
223 * the trajectory and we have already read the line in the *
224 * previous call (VERY DIRTY). */
225 bFinished = FALSE;
226 do {
227 bTime = (strcmp(line,"TIMESTEP") == 0);
228 bAtoms = (strcmp(line,"POSITION") == 0);
229 bPos = (bAtoms || (strcmp(line,"POSITIONRED") == 0));
230 bVel = (strncmp(line,"VELOCITY",8) == 0);
231 bBox = (strcmp(line,"BOX") == 0);
232 if (bTime) {
233 if (!fr->bTime && !fr->bX) {
234 fr->bStep = bTime;
235 fr->bTime = bTime;
236 do
237 bFinished = (fgets2(line,STRLEN,fp) == NULL);
238 while (!bFinished && (line[0] == '#'));
239 sscanf(line,"%15d%15lf",&(fr->step),&db1);
240 fr->time = db1;
241 } else
242 bFinished = TRUE;
243 }
244 if (bPos) {
245 if (!fr->bX) {
246 fr->bAtoms = bAtoms;
247 fr->bX = bPos;
248 natoms = read_g96_pos(line,symtab,fp,infile,fr);
249 } else
250 bFinished = TRUE;
251 }
252 if (fr->v && bVel) {
253 fr->bV = bVel;
254 natoms = read_g96_vel(line,fp,infile,fr);
255 }
256 if (bBox) {
257 fr->bBox = bBox;
258 clear_mat(fr->box);
259 bEnd = FALSE;
260 while (!bEnd && fgets2(line,STRLEN,fp)) {
261 bEnd = (strncmp(line,"END",3) == 0);
262 if (!bEnd && (line[0] != '#')) {
263 nbp = sscanf(line,"%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf",
264 &db1,&db2,&db3,&db4,&db5,&db6,&db7,&db8,&db9);
265 if (nbp < 3)
266 gmx_fatal(FARGS,"Found a BOX line, but no box in %s",infile);
267 fr->box[XX][XX] = db1;
268 fr->box[YY][YY] = db2;
269 fr->box[ZZ][ZZ] = db3;
270 if (nbp == 9) {
271 fr->box[XX][YY] = db4;
272 fr->box[XX][ZZ] = db5;
273 fr->box[YY][XX] = db6;
274 fr->box[YY][ZZ] = db7;
275 fr->box[ZZ][XX] = db8;
276 fr->box[ZZ][YY] = db9;
277 }
278 }
279 }
280 bFinished = TRUE;
281 }
282 } while (!bFinished && fgets2(line,STRLEN,fp));
283
284 free_symtab(symtab);
285
286 fr->natoms = natoms;
287
288 return natoms;
289 }
290
291 void write_g96_conf(FILE *out,t_trxframe *fr,
292 int nindex,atom_id *index)
293 {
294 t_atoms *atoms;
295 int nout,i,a;
296
297 atoms = fr->atoms;
298
299 if (index)
300 nout = nindex;
301 else
302 nout = fr->natoms;
303
304 if (fr->bTitle)
305 fprintf(out,"TITLE\n%s\nEND\n",fr->title);
306 if (fr->bStep || fr->bTime)
307 /* Officially the time format is %15.9, which is not enough for 10 ns */
308 fprintf(out,"TIMESTEP\n%15d%15.6f\nEND\n",fr->step,fr->time);
309 if (fr->bX) {
310 if (fr->bAtoms) {
311 fprintf(out,"POSITION\n");
312 for(i=0; i<nout; i++) {
313 if (index) a = index[i]; else a = i;
314 fprintf(out,"%5d %-5s %-5s%7d%15.9f%15.9f%15.9f\n",
315 (atoms->resinfo[atoms->atom[a].resind].nr) % 100000,
316 *atoms->resinfo[atoms->atom[a].resind].name,
317 *atoms->atomname[a],(i+1) % 10000000,
318 fr->x[a][XX],fr->x[a][YY],fr->x[a][ZZ]);
319 }
320 } else {
321 fprintf(out,"POSITIONRED\n");
322 for(i=0; i<nout; i++) {
323 if (index) a = index[i]; else a = i;
324 fprintf(out,"%15.9f%15.9f%15.9f\n",
325 fr->x[a][XX],fr->x[a][YY],fr->x[a][ZZ]);
326 }
327 }
328 fprintf(out,"END\n");
329 }
330 if (fr->bV) {
331 if (fr->bAtoms) {
332 fprintf(out,"VELOCITY\n");
333 for(i=0; i<nout; i++) {
334 if (index) a = index[i]; else a = i;
335 fprintf(out,"%5d %-5s %-5s%7d%15.9f%15.9f%15.9f\n",
336 (atoms->resinfo[atoms->atom[a].resind].nr) % 100000,
337 *atoms->resinfo[atoms->atom[a].resind].name,
338 *atoms->atomname[a],(i+1) % 10000000,
339 fr->v[a][XX],fr->v[a][YY],fr->v[a][ZZ]);
340 }
341 } else {
342 fprintf(out,"VELOCITYRED\n");
343 for(i=0; i<nout; i++) {
344 if (index) a = index[i]; else a = i;
345 fprintf(out,"%15.9f%15.9f%15.9f\n",
346 fr->v[a][XX],fr->v[a][YY],fr->v[a][ZZ]);
347 }
348 }
349 fprintf(out,"END\n");
350 }
351 if (fr->bBox) {
352 fprintf(out,"BOX\n");
353 fprintf(out,"%15.9f%15.9f%15.9f",
354 fr->box[XX][XX],fr->box[YY][YY],fr->box[ZZ][ZZ]);
355 if (fr->box[XX][YY] || fr->box[XX][ZZ] || fr->box[YY][XX] ||
356 fr->box[YY][ZZ] || fr->box[ZZ][XX] ||fr->box[ZZ][YY])
357 fprintf(out,"%15.9f%15.9f%15.9f%15.9f%15.9f%15.9f",
358 fr->box[XX][YY],fr->box[XX][ZZ],fr->box[YY][XX],
359 fr->box[YY][ZZ],fr->box[ZZ][XX],fr->box[ZZ][YY]);
360 fprintf(out,"\n");
361 fprintf(out,"END\n");
362 }
363 }
364
365 static int get_espresso_word(FILE *fp,char word[])
366 {
367 int ret,nc,i;
368
369 ret = 0;
370 nc = 0;
371
372 do {
373 i = fgetc(fp);
374 if (i != EOF) {
375 if (i == ' ' || i == '\n' || i == '\t') {
376 if (nc > 0)
377 ret = 1;
378 } else if (i == '{') {
379 if (nc == 0)
380 word[nc++] = '{';
381 ret = 2;
382 } else if (i == '}') {
383 if (nc == 0)
384 word[nc++] = '}';
385 ret = 3;
386 } else {
387 word[nc++] = (char)i;
388 }
389 }
390 } while (i != EOF && ret == 0);
391
392 word[nc] = '\0';
393
394 /* printf("'%s'\n",word); */
395
396 return ret;
397 }
398
399 static int check_open_parenthesis(FILE *fp,int r,
400 const char *infile,const char *keyword)
401 {
402 int level_inc;
403 char word[STRLEN];
404
405 level_inc = 0;
406 if (r == 2) {
407 level_inc++;
408 } else {
409 r = get_espresso_word(fp,word);
410 if (r == 2)
411 level_inc++;
412 else
413 gmx_fatal(FARGS,"Expected '{' after '%s' in file '%s'",
414 keyword,infile);
415 }
416
417 return level_inc;
418 }
419
420 static int check_close_parenthesis(FILE *fp,int r,
421 const char *infile,const char *keyword)
422 {
423 int level_inc;
424 char word[STRLEN];
425
426 level_inc = 0;
427 if (r == 3) {
428 level_inc--;
429 } else {
430 r = get_espresso_word(fp,word);
431 if (r == 3)
432 level_inc--;
433 else
434 gmx_fatal(FARGS,"Expected '}' after section '%s' in file '%s'",
435 keyword,infile);
436 }
437
438 return level_inc;
439 }
440
441 enum { espID, espPOS, espTYPE, espQ, espV, espF, espMOLECULE, espNR };
442 const char *esp_prop[espNR] = { "id", "pos", "type", "q", "v", "f",
443 "molecule" };
444
445 static void read_espresso_conf(const char *infile,
446 t_atoms *atoms,rvec x[],rvec *v,matrix box)
447 {
448 t_symtab *symtab=NULL;
449 FILE *fp;
450 char word[STRLEN],buf[STRLEN];
451 int natoms,level,npar,r,nprop,p,i,m,molnr;
452 int prop[32];
453 double d;
454 gmx_bool bFoundParticles,bFoundProp,bFoundVariable,bMol;
455
456 if (!symtab) {
457 snew(symtab,1);
458 open_symtab(symtab);
459 }
460
461 clear_mat(box);
462
463 fp = gmx_fio_fopen(infile,"r");
464
465 bFoundParticles = FALSE;
466 bFoundVariable = FALSE;
467 bMol = FALSE;
468 level = 0;
469 while ((r=get_espresso_word(fp,word))) {
470 if (level==1 && strcmp(word,"particles")==0 && !bFoundParticles) {
471 bFoundParticles = TRUE;
472 level += check_open_parenthesis(fp,r,infile,"particles");
473 nprop = 0;
474 while (level == 2 && (r=get_espresso_word(fp,word))) {
475 bFoundProp = FALSE;
476 for(p=0; p<espNR; p++) {
477 if (strcmp(word,esp_prop[p]) == 0) {
478 bFoundProp = TRUE;
479 prop[nprop++] = p;
480 /* printf(" prop[%d] = %s\n",nprop-1,esp_prop[prop[nprop-1]]); */
481 }
482 }
483 if (!bFoundProp && word[0] != '}') {
484 gmx_fatal(FARGS,"Can not read Espresso files with particle property '%s'",word);
485 }
486 if (bFoundProp && p == espMOLECULE)
487 bMol = TRUE;
488 if (r == 3)
489 level--;
490 }
491
492 i = 0;
493 while (level > 0 && (r=get_espresso_word(fp,word))) {
494 if (r == 2) {
495 level++;
496 } else if (r == 3) {
497 level--;
498 }
499 if (level == 2) {
500 for(p=0; p<nprop; p++) {
501 switch (prop[p]) {
502 case espID:
503 r = get_espresso_word(fp,word);
504 /* Not used */
505 break;
506 case espPOS:
507 for(m=0; m<3; m++) {
508 r = get_espresso_word(fp,word);
509 sscanf(word,"%lf",&d);
510 x[i][m] = d;
511 }
512 break;
513 case espTYPE:
514 r = get_espresso_word(fp,word);
515 atoms->atom[i].type = strtol(word, NULL, 10);
516 break;
517 case espQ:
518 r = get_espresso_word(fp,word);
519 sscanf(word,"%lf",&d);
520 atoms->atom[i].q = d;
521 break;
522 case espV:
523 for(m=0; m<3; m++) {
524 r = get_espresso_word(fp,word);
525 sscanf(word,"%lf",&d);
526 v[i][m] = d;
527 }
528 break;
529 case espF:
530 for(m=0; m<3; m++) {
531 r = get_espresso_word(fp,word);
532 /* not used */
533 }
534 break;
535 case espMOLECULE:
536 r = get_espresso_word(fp,word);
537 molnr = strtol(word, NULL, 10);
538 if (i == 0 ||
539 atoms->resinfo[atoms->atom[i-1].resind].nr != molnr) {
540 atoms->atom[i].resind =
541 (i == 0 ? 0 : atoms->atom[i-1].resind+1);
542 atoms->resinfo[atoms->atom[i].resind].nr = molnr;
543 atoms->resinfo[atoms->atom[i].resind].ic = ' ';
544 atoms->resinfo[atoms->atom[i].resind].chainid = ' ';
545 atoms->resinfo[atoms->atom[i].resind].chainnum = molnr; /* Not sure if this is right? */
546 } else {
547 atoms->atom[i].resind = atoms->atom[i-1].resind;
548 }
549 break;
550 }
551 }
552 /* Generate an atom name from the particle type */
553 sprintf(buf,"T%d",atoms->atom[i].type);
554 atoms->atomname[i] = put_symtab(symtab,buf);
555 if (bMol) {
556 if (i == 0 || atoms->atom[i].resind != atoms->atom[i-1].resind) {
557 atoms->resinfo[atoms->atom[i].resind].name =
558 put_symtab(symtab,"MOL");
559 }
560 } else {
561 /* Residue number is the atom number */
562 atoms->atom[i].resind = i;
563 /* Generate an residue name from the particle type */
564 if (atoms->atom[i].type < 26) {
565 sprintf(buf,"T%c",'A'+atoms->atom[i].type);
566 } else {
567 sprintf(buf,"T%c%c",
568 'A'+atoms->atom[i].type/26,'A'+atoms->atom[i].type%26);
569 }
570 t_atoms_set_resinfo(atoms,i,symtab,buf,i,' ',0,' ');
571 }
572
573 if (r == 3)
574 level--;
575 i++;
576 }
577 }
578 atoms->nres = atoms->nr;
579
580 if (i != atoms->nr) {
581 gmx_fatal(FARGS,"Internal inconsistency in Espresso routines, read %d atoms, expected %d atoms",i,atoms->nr);
582 }
583 } else if (level==1 && strcmp(word,"variable")==0 && !bFoundVariable) {
584 bFoundVariable = TRUE;
585 level += check_open_parenthesis(fp,r,infile,"variable");
586 while (level==2 && (r=get_espresso_word(fp,word))) {
587 if (level==2 && strcmp(word,"box_l") == 0) {
588 for(m=0; m<3; m++) {
589 r = get_espresso_word(fp,word);
590 sscanf(word,"%lf",&d);
591 box[m][m] = d;
592 }
593 level += check_close_parenthesis(fp,r,infile,"box_l");
594 }
595 }
596 } else if (r == 2) {
597 level++;
598 } else if (r == 3) {
599 level--;
600 }
601 }
602
603 if (!bFoundParticles) {
604 fprintf(stderr,"Did not find a particles section in Espresso file '%s'\n",
605 infile);
606 }
607
608 gmx_fio_fclose(fp);
609 }
610
611 static int get_espresso_coordnum(const char *infile)
612 {
613 FILE *fp;
614 char word[STRLEN];
615 int natoms,level,r;
616 gmx_bool bFoundParticles;
617
618 natoms = 0;
619
620 fp = gmx_fio_fopen(infile,"r");
621
622 bFoundParticles = FALSE;
623 level = 0;
624 while ((r=get_espresso_word(fp,word)) && !bFoundParticles) {
625 if (level==1 && strcmp(word,"particles")==0 && !bFoundParticles) {
626 bFoundParticles = TRUE;
627 level += check_open_parenthesis(fp,r,infile,"particles");
628 while (level > 0 && (r=get_espresso_word(fp,word))) {
629 if (r == 2) {
630 level++;
631 if (level == 2)
632 natoms++;
633 } else if (r == 3) {
634 level--;
635 }
636 }
637 } else if (r == 2) {
638 level++;
639 } else if (r == 3) {
640 level--;
641 }
642 }
643 if (!bFoundParticles) {
644 fprintf(stderr,"Did not find a particles section in Espresso file '%s'\n",
645 infile);
646 }
647
648 gmx_fio_fclose(fp);
649
650 return natoms;
651 }
652
653 static void write_espresso_conf_indexed(FILE *out,const char *title,
654 t_atoms *atoms,int nx,atom_id *index,
655 rvec *x,rvec *v,matrix box)
656 {
657 int i,j;
658
659 fprintf(out,"# %s\n",title);
660 if (TRICLINIC(box)) {
661 gmx_warning("The Espresso format does not support triclinic unit-cells");
662 }
663 fprintf(out,"{variable {box_l %f %f %f}}\n",box[0][0],box[1][1],box[2][2]);
664
665 fprintf(out,"{particles {id pos type q%s}\n",v ? " v" : "");
666 for(i=0; i<nx; i++) {
667 if (index)
668 j = index[i];
669 else
670 j = i;
671 fprintf(out,"\t{%d %f %f %f %d %g",
672 j,x[j][XX],x[j][YY],x[j][ZZ],
673 atoms->atom[j].type,atoms->atom[j].q);
674 if (v)
675 fprintf(out," %f %f %f",v[j][XX],v[j][YY],v[j][ZZ]);
676 fprintf(out,"}\n");
677 }
678 fprintf(out,"}\n");
679 }
680
681 static void get_coordnum_fp (FILE *in, char *title, int *natoms)
682 {
683 char line[STRLEN+1];
684
685 fgets2 (title,STRLEN,in);
686 fgets2 (line,STRLEN,in);
687 if (sscanf (line,"%d",natoms) != 1) {
688 gmx_fatal(FARGS,"gro file does not have the number of atoms on the second line");
689 }
690 }
691
692 static void get_coordnum (const char *infile,int *natoms)
693 {
694 FILE *in;
695 char title[STRLEN];
696
697 in=gmx_fio_fopen(infile,"r");
698 get_coordnum_fp(in,title,natoms);
699 gmx_fio_fclose (in);
700 }
701
702 static gmx_bool get_w_conf(FILE *in,const char *infile,char *title,
703 t_atoms *atoms, int *ndec, rvec x[],rvec *v, matrix box)
704 {
705 t_symtab *symtab=NULL;
706 char name[6];
707 char line[STRLEN+1],*ptr;
708 char buf[256];
709 double x1,y1,z1,x2,y2,z2;
710 rvec xmin,xmax;
711 int natoms,i,m,resnr,newres,oldres,ddist,c;
712 gmx_bool bFirst,bVel;
713 char *p1,*p2,*p3;
714
715 newres = -1;
716 oldres = NOTSET; /* Unlikely number for the first residue! */
717 ddist = 0;
718
719 if (!symtab) {
720 snew(symtab,1);
721 open_symtab(symtab);
722 }
723
724 /* Read the title and number of atoms */
725 get_coordnum_fp(in,title,&natoms);
726
727 if (natoms > atoms->nr)
728 gmx_fatal(FARGS,"gro file contains more atoms (%d) than expected (%d)",
729 natoms,atoms->nr);
730 else if (natoms < atoms->nr)
731 fprintf(stderr,"Warning: gro file contains less atoms (%d) than expected"
732 " (%d)\n",natoms,atoms->nr);
733
734 bFirst=TRUE;
735
736 bVel = FALSE;
737
738 /* just pray the arrays are big enough */
739 for (i=0; (i < natoms) ; i++) {
740 if ((fgets2 (line,STRLEN,in)) == NULL) {
741 unexpected_eof(infile,i+2);
742 }
743 if (strlen(line) < 39)
744 gmx_fatal(FARGS,"Invalid line in %s for atom %d:\n%s",infile,i+1,line);
745
746 /* determine read precision from distance between periods
747 (decimal points) */
748 if (bFirst) {
749 bFirst=FALSE;
750 p1=strchr(line,'.');
751 if (p1 == NULL)
752 gmx_fatal(FARGS,"A coordinate in file %s does not contain a '.'",infile);
753 p2=strchr(&p1[1],'.');
754 if (p2 == NULL)
755 gmx_fatal(FARGS,"A coordinate in file %s does not contain a '.'",infile);
756 ddist = p2 - p1;
757 *ndec = ddist - 5;
758
759 p3=strchr(&p2[1],'.');
760 if (p3 == NULL)
761 gmx_fatal(FARGS,"A coordinate in file %s does not contain a '.'",infile);
762
763 if (p3 - p2 != ddist)
764 gmx_fatal(FARGS,"The spacing of the decimal points in file %s is not consistent for x, y and z",infile);
765 }
766
767 /* residue number*/
768 memcpy(name,line,5);
769 name[5]='\0';
770 sscanf(name,"%d",&resnr);
771 memcpy(name,line+5,5);
772 name[5]='\0';
773 if (resnr != oldres) {
774 oldres = resnr;
775 newres++;
776 if (newres >= natoms)
777 gmx_fatal(FARGS,"More residues than atoms in %s (natoms = %d)",
778 infile,natoms);
779 atoms->atom[i].resind = newres;
780 t_atoms_set_resinfo(atoms,i,symtab,name,resnr,' ',0,' ');
781 } else {
782 atoms->atom[i].resind = newres;
783 }
784
785 /* atomname */
786 memcpy(name,line+10,5);
787 atoms->atomname[i]=put_symtab(symtab,name);
788
789 /* eventueel controle atomnumber met i+1 */
790
791 /* coordinates (start after residue shit) */
792 ptr = line + 20;
793 /* Read fixed format */
794 for(m=0; m<DIM; m++) {
795 for(c=0; (c<ddist && ptr[0]); c++) {
796 buf[c] = ptr[0];
797 ptr++;
798 }
799 buf[c] = '\0';
800 if (sscanf (buf,"%lf %lf",&x1,&x2) != 1) {
801 gmx_fatal(FARGS,"Something is wrong in the coordinate formatting of file %s. Note that gro is fixed format (see the manual)",infile);
802 } else {
803 x[i][m] = x1;
804 }
805 }
806
807 /* velocities (start after residues and coordinates) */
808 if (v) {
809 /* Read fixed format */
810 for(m=0; m<DIM; m++) {
811 for(c=0; (c<ddist && ptr[0]); c++) {
812 buf[c] = ptr[0];
813 ptr++;
814 }
815 buf[c] = '\0';
816 if (sscanf (buf,"%lf",&x1) != 1) {
817 v[i][m] = 0;
818 } else {
819 v[i][m] = x1;
820 bVel = TRUE;
821 }
822 }
823 }
824 }
825 atoms->nres = newres + 1;
826
827 /* box */
828 fgets2 (line,STRLEN,in);
829 if (sscanf (line,"%lf%lf%lf",&x1,&y1,&z1) != 3) {
830 gmx_warning("Bad box in file %s",infile);
831
832 /* Generate a cubic box */
833 for(m=0; (m<DIM); m++)
834 xmin[m]=xmax[m]=x[0][m];
835 for(i=1; (i<atoms->nr); i++)
836 for(m=0; (m<DIM); m++) {
837 xmin[m]=min(xmin[m],x[i][m]);
838 xmax[m]=max(xmax[m],x[i][m]);
839 }
840 for (i=0; i<DIM; i++) for (m=0; m<DIM; m++) box[i][m]=0.0;
841 for(m=0; (m<DIM); m++)
842 box[m][m]=(xmax[m]-xmin[m]);
843 fprintf(stderr,"Generated a cubic box %8.3f x %8.3f x %8.3f\n",
844 box[XX][XX],box[YY][YY],box[ZZ][ZZ]);
845 }
846 else {
847 /* We found the first three values, the diagonal elements */
848 box[XX][XX]=x1;
849 box[YY][YY]=y1;
850 box[ZZ][ZZ]=z1;
851 if (sscanf (line,"%*f%*f%*f%lf%lf%lf%lf%lf%lf",
852 &x1,&y1,&z1,&x2,&y2,&z2) != 6)
853 x1=y1=z1=x2=y2=z2=0.0;
854 box[XX][YY] = x1;
855 box[XX][ZZ] = y1;
856 box[YY][XX] = z1;
857 box[YY][ZZ] = x2;
858 box[ZZ][XX] = y2;
859 box[ZZ][YY] = z2;
860 }
861 close_symtab(symtab);
862
863 return bVel;
864 }
865
866 static void read_whole_conf(const char *infile,char *title,
867 t_atoms *atoms, rvec x[],rvec *v, matrix box)
868 {
869 FILE *in;
870 int ndec;
871
872 /* open file */
873 in=gmx_fio_fopen(infile,"r");
874
875 get_w_conf(in, infile, title, atoms, &ndec, x, v, box);
876
877 gmx_fio_fclose(in);
878 }
879
880 static void get_conf(FILE *in,char *title,int *natoms,
881 rvec x[],rvec *v,matrix box)
882 {
883 t_atoms atoms;
884 int ndec;
885
886 atoms.nr=*natoms;
887 snew(atoms.atom,*natoms);
888 atoms.nres=*natoms;
889 snew(atoms.resinfo,*natoms);
890 snew(atoms.atomname,*natoms);
891
892 get_w_conf(in,title,title,&atoms,&ndec,x,v,box);
893
894 sfree(atoms.atom);
895 sfree(atoms.resinfo);
896 sfree(atoms.atomname);
897 }
898
899 gmx_bool gro_next_x_or_v(FILE *status,t_trxframe *fr)
900 {
901 t_atoms atoms;
902 char title[STRLEN],*p;
903 double tt;
904 int ndec,i;
905
906 if (gmx_eof(status))
907 return FALSE;
908
909 atoms.nr=fr->natoms;
910 snew(atoms.atom,fr->natoms);
911 atoms.nres=fr->natoms;
912 snew(atoms.resinfo,fr->natoms);
913 snew(atoms.atomname,fr->natoms);
914
915 fr->bV = get_w_conf(status,title,title,&atoms,&ndec,fr->x,fr->v,fr->box);
916 fr->bPrec = TRUE;
917 fr->prec = 1;
918 /* prec = 10^ndec: */
919 for(i=0; i<ndec; i++)
920 fr->prec *= 10;
921 fr->title = title;
922 fr->bTitle = TRUE;
923 fr->bX = TRUE;
924 fr->bBox = TRUE;
925
926 sfree(atoms.atom);
927 sfree(atoms.resinfo);
928 sfree(atoms.atomname);
929
930 if ((p=strstr(title,"t=")) != NULL) {
931 p+=2;
932 if (sscanf(p,"%lf",&tt)==1) {
933 fr->time = tt;
934 fr->bTime = TRUE;
935 } else {
936 fr->time = 0;
937 fr->bTime = FALSE;
938 }
939 }
940
941 if (atoms.nr != fr->natoms)
942 gmx_fatal(FARGS,"Number of atoms in gro frame (%d) doesn't match the number in the previous frame (%d)",atoms.nr,fr->natoms);
943
944 return TRUE;
945 }
946
947 int gro_first_x_or_v(FILE *status,t_trxframe *fr)
948 {
949 char title[STRLEN];
950
951 frewind(status);
952 fprintf(stderr,"Reading frames from gro file");
953 get_coordnum_fp(status, title, &fr->natoms);
954 frewind(status);
955 fprintf(stderr," '%s', %d atoms.\n",title, fr->natoms);
956 fr->bTitle = TRUE;
957 fr->title = title;
958 if (fr->natoms==0)
959 gmx_file("No coordinates in gro file");
960
961 snew(fr->x,fr->natoms);
962 snew(fr->v,fr->natoms);
963 gro_next_x_or_v(status, fr);
964
965 return fr->natoms;
966 }
967
968 static void make_hconf_format(int pr,gmx_bool bVel,char format[])
969 {
970 int l,vpr;
971
972 /* build format string for printing,
973 something like "%8.3f" for x and "%8.4f" for v */
974 if (pr<0)
975 pr=0;
976 if (pr>30)
977 pr=30;
978 l=pr+5;
979 vpr=pr+1;
980 if (bVel)
981 sprintf(format,"%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df\n",
982 l,pr,l,pr,l,pr,l,vpr,l,vpr,l,vpr);
983 else
984 sprintf(format,"%%%d.%df%%%d.%df%%%d.%df\n",l,pr,l,pr,l,pr);
985
986 }
987
988 static void write_hconf_box(FILE *out,int pr,matrix box)
989 {
990 char format[100];
991 int l;
992
993 if (pr<5)
994 pr=5;
995 l=pr+5;
996
997 if (box[XX][YY] || box[XX][ZZ] || box[YY][XX] || box[YY][ZZ] ||
998 box[ZZ][XX] || box[ZZ][YY]) {
999 sprintf(format,"%%%d.%df%%%d.%df%%%d.%df"
1000 "%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df\n",
1001 l,pr,l,pr,l,pr,l,pr,l,pr,l,pr,l,pr,l,pr,l,pr);
1002 fprintf(out,format,
1003 box[XX][XX],box[YY][YY],box[ZZ][ZZ],
1004 box[XX][YY],box[XX][ZZ],box[YY][XX],
1005 box[YY][ZZ],box[ZZ][XX],box[ZZ][YY]);
1006 } else {
1007 sprintf(format,"%%%d.%df%%%d.%df%%%d.%df\n",l,pr,l,pr,l,pr);
1008 fprintf(out,format,
1009 box[XX][XX],box[YY][YY],box[ZZ][ZZ]);
1010 }
1011 }
1012
1013 void write_hconf_indexed_p(FILE *out,const char *title,t_atoms *atoms,
1014 int nx,atom_id index[], int pr,
1015 rvec *x,rvec *v,matrix box)
1016 {
1017 char resnm[6],nm[6],format[100];
1018 int ai,i,resind,resnr;
1019
1020 bromacs(format,99);
1021 fprintf (out,"%s\n",(title && title[0])?title:format);
1022 fprintf (out,"%5d\n",nx);
1023
1024 make_hconf_format(pr,v!=NULL,format);
1025
1026 for (i=0; (i<nx); i++) {
1027 ai=index[i];
1028
1029 resind = atoms->atom[ai].resind;
1030 strncpy(resnm," ??? ",sizeof(resnm)-1);
1031 if (resind < atoms->nres) {
1032 strncpy(resnm,*atoms->resinfo[resind].name,sizeof(resnm)-1);
1033 resnr = atoms->resinfo[resind].nr;
1034 } else {
1035 strncpy(resnm," ??? ",sizeof(resnm)-1);
1036 resnr = resind + 1;
1037 }
1038
1039 if (atoms->atom)
1040 strncpy(nm,*atoms->atomname[ai],sizeof(nm)-1);
1041 else
1042 strncpy(nm," ??? ",sizeof(nm)-1);
1043
1044 fprintf(out,"%5d%-5.5s%5.5s%5d",resnr%100000,resnm,nm,(ai+1)%100000);
1045 /* next fprintf uses built format string */
1046 if (v)
1047 fprintf(out,format,
1048 x[ai][XX], x[ai][YY], x[ai][ZZ], v[ai][XX],v[ai][YY],v[ai][ZZ]);
1049 else
1050 fprintf(out,format,
1051 x[ai][XX], x[ai][YY], x[ai][ZZ]);
1052 }
1053
1054 write_hconf_box(out,pr,box);
1055
1056 fflush(out);
1057 }
1058
1059 static void write_hconf_mtop(FILE *out,const char *title,gmx_mtop_t *mtop,
1060 int pr,
1061 rvec *x,rvec *v,matrix box)
1062 {
1063 char format[100];
1064 int i,resnr;
1065 gmx_mtop_atomloop_all_t aloop;
1066 t_atom *atom;
1067 char *atomname,*resname;
1068
1069 bromacs(format,99);
1070 fprintf (out,"%s\n",(title && title[0])?title:format);
1071 fprintf (out,"%5d\n",mtop->natoms);
1072
1073 make_hconf_format(pr,v!=NULL,format);
1074
1075 aloop = gmx_mtop_atomloop_all_init(mtop);
1076 while (gmx_mtop_atomloop_all_next(aloop,&i,&atom)) {
1077 gmx_mtop_atomloop_all_names(aloop,&atomname,&resnr,&resname);
1078
1079 fprintf(out,"%5d%-5.5s%5.5s%5d",
1080 resnr%100000,resname,atomname,(i+1)%100000);
1081 /* next fprintf uses built format string */
1082 if (v)
1083 fprintf(out,format,
1084 x[i][XX], x[i][YY], x[i][ZZ], v[i][XX],v[i][YY],v[i][ZZ]);
1085 else
1086 fprintf(out,format,
1087 x[i][XX], x[i][YY], x[i][ZZ]);
1088 }
1089
1090 write_hconf_box(out,pr,box);
1091
1092 fflush(out);
1093 }
1094
1095 void write_hconf_p(FILE *out,const char *title,t_atoms *atoms, int pr,
1096 rvec *x,rvec *v,matrix box)
1097 {
1098 atom_id *aa;
1099 int i;
1100
1101 snew(aa,atoms->nr);
1102 for(i=0; (i<atoms->nr); i++)
1103 aa[i]=i;
1104 write_hconf_indexed_p(out,title,atoms,atoms->nr,aa,pr,x,v,box);
1105 sfree(aa);
1106 }
1107
1108 void write_conf_p(const char *outfile, const char *title,
1109 t_atoms *atoms, int pr,
1110 rvec *x, rvec *v,matrix box)
1111 {
1112 FILE *out;
1113
1114 out=gmx_fio_fopen(outfile,"w");
1115 write_hconf_p(out,title,atoms,pr,x,v,box);
1116
1117 gmx_fio_fclose (out);
1118 }
1119
1120 static void write_conf(const char *outfile, const char *title, t_atoms *atoms,
1121 rvec *x, rvec *v,matrix box)
1122 {
1123 write_conf_p(outfile, title, atoms, 3, x, v, box);
1124 }
1125
1126 void write_sto_conf_indexed(const char *outfile,const char *title,
1127 t_atoms *atoms,
1128 rvec x[],rvec *v,int ePBC,matrix box,
1129 atom_id nindex,atom_id index[])
1130 {
1131 FILE *out;
1132 int ftp;
1133 t_trxframe fr;
1134
1135 ftp=fn2ftp(outfile);
1136 switch (ftp) {
1137 case efGRO:
1138 out=gmx_fio_fopen(outfile,"w");
1139 write_hconf_indexed_p(out, title, atoms, nindex, index, 3, x, v, box);
1140 gmx_fio_fclose(out);
1141 break;
1142 case efG96:
1143 clear_trxframe(&fr,TRUE);
1144 fr.bTitle = TRUE;
1145 fr.title = title;
1146 fr.natoms = atoms->nr;
1147 fr.bAtoms = TRUE;
1148 fr.atoms = atoms;
1149 fr.bX = TRUE;
1150 fr.x = x;
1151 if (v) {
1152 fr.bV = TRUE;
1153 fr.v = v;
1154 }
1155 fr.bBox = TRUE;
1156 copy_mat(box,fr.box);
1157 out=gmx_fio_fopen(outfile,"w");
1158 write_g96_conf(out, &fr, nindex, index);
1159 gmx_fio_fclose(out);
1160 break;
1161 case efPDB:
1162 case efBRK:
1163 case efENT:
1164 case efPQR:
1165 out=gmx_fio_fopen(outfile,"w");
1166 write_pdbfile_indexed(out,title,atoms,x,ePBC,box,' ',-1,nindex,index,NULL,TRUE);
1167 gmx_fio_fclose(out);
1168 break;
1169 case efESP:
1170 out=gmx_fio_fopen(outfile,"w");
1171 write_espresso_conf_indexed(out, title, atoms, nindex, index, x, v, box);
1172 gmx_fio_fclose(out);
1173 break;
1174 case efTPR:
1175 case efTPB:
1176 case efTPA:
1177 gmx_fatal(FARGS,"Sorry, can not write a topology to %s",outfile);
1178 break;
1179 default:
1180 gmx_incons("Not supported in write_sto_conf_indexed");
1181 }
1182 }
1183
1184 static void write_xyz_conf(const char *outfile,const char *title,
1185 t_atoms *atoms,rvec *x)
1186 {
1187 FILE *fp;
1188 int i,anr;
1189 real value;
1190 char *ptr,*name;
1191 gmx_atomprop_t aps=gmx_atomprop_init();
1192
1193 fp = gmx_fio_fopen(outfile,"w");
1194 fprintf(fp,"%3d\n",atoms->nr);
1195 fprintf(fp,"%s\n",title);
1196 for(i=0; (i<atoms->nr); i++) {
1197 anr = atoms->atom[i].atomnumber;
1198 name = *atoms->atomname[i];
1199 if (anr == NOTSET) {
1200 if (gmx_atomprop_query(aps,epropElement,"???",name,&value))
1201 anr = gmx_nint(value);
1202 }
1203 if ((ptr = gmx_atomprop_element(aps,anr)) == NULL)
1204 ptr = name;
1205 fprintf(fp,"%3s%10.5f%10.5f%10.5f\n",ptr,
1206 10*x[i][XX],10*x[i][YY],10*x[i][ZZ]);
1207 }
1208 gmx_fio_fclose(fp);
1209 gmx_atomprop_destroy(aps);
1210 }
1211
1212 void write_sto_conf(const char *outfile,const char *title,t_atoms *atoms,
1213 rvec x[],rvec *v,int ePBC,matrix box)
1214 {
1215 FILE *out;
1216 int ftp;
1217 t_trxframe fr;
1218
1219 ftp=fn2ftp(outfile);
1220 switch (ftp) {
1221 case efGRO:
1222 write_conf(outfile, title, atoms, x, v, box);
1223 break;
1224 case efG96:
1225 clear_trxframe(&fr,TRUE);
1226 fr.bTitle = TRUE;
1227 fr.title = title;
1228 fr.natoms = atoms->nr;
1229 fr.bAtoms = TRUE;
1230 fr.atoms = atoms;
1231 fr.bX = TRUE;
1232 fr.x = x;
1233 if (v) {
1234 fr.bV = TRUE;
1235 fr.v = v;
1236 }
1237 fr.bBox = TRUE;
1238 copy_mat(box,fr.box);
1239 out=gmx_fio_fopen(outfile,"w");
1240 write_g96_conf(out, &fr, -1, NULL);
1241 gmx_fio_fclose(out);
1242 break;
1243 case efXYZ:
1244 write_xyz_conf(outfile,(strlen(title) > 0) ? title : outfile,atoms,x);
1245 break;
1246 case efPDB:
1247 case efBRK:
1248 case efENT:
1249 out=gmx_fio_fopen(outfile,"w");
1250 write_pdbfile(out, title, atoms, x, ePBC, box, ' ', -1,NULL,TRUE);
1251 gmx_fio_fclose(out);
1252 break;
1253 case efESP:
1254 out=gmx_fio_fopen(outfile,"w");
1255 write_espresso_conf_indexed(out, title, atoms, atoms->nr, NULL, x, v, box);
1256 gmx_fio_fclose(out);
1257 break;
1258 case efTPR:
1259 case efTPB:
1260 case efTPA:
1261 gmx_fatal(FARGS,"Sorry, can not write a topology to %s",outfile);
1262 break;
1263 default:
1264 gmx_incons("Not supported in write_sto_conf");
1265 }
1266 }
1267
1268 void write_sto_conf_mtop(const char *outfile,const char *title,
1269 gmx_mtop_t *mtop,
1270 rvec x[],rvec *v,int ePBC,matrix box)
1271 {
1272 int ftp;
1273 FILE *out;
1274 t_atoms atoms;
1275
1276 ftp=fn2ftp(outfile);
1277 switch (ftp) {
1278 case efGRO:
1279 out = gmx_fio_fopen(outfile,"w");
1280 write_hconf_mtop(out,title,mtop,3,x,v,box);
1281 gmx_fio_fclose(out);
1282 break;
1283 default:
1284 /* This is a brute force approach which requires a lot of memory.
1285 * We should implement mtop versions of all writing routines.
1286 */
1287 atoms = gmx_mtop_global_atoms(mtop);
1288
1289 write_sto_conf(outfile,title,&atoms,x,v,ePBC,box);
1290
1291 done_atom(&atoms);
1292 break;
1293 }
1294 }
1295
1296 static int get_xyz_coordnum(const char *infile)
1297 {
1298 FILE *fp;
1299 int n;
1300
1301 fp = gmx_fio_fopen(infile,"r");
1302 if (fscanf(fp,"%d",&n) != 1)
1303 gmx_fatal(FARGS,"Can not read number of atoms from %s",infile);
1304 gmx_fio_fclose(fp);
1305
1306 return n;
1307 }
1308
1309 static void read_xyz_conf(const char *infile,char *title,
1310 t_atoms *atoms,rvec *x)
1311 {
1312 FILE *fp;
1313 int i,n;
1314 double xx,yy,zz;
1315 t_symtab *tab;
1316 char atomnm[32],buf[STRLEN];
1317
1318 snew(tab,1);
1319 fp = gmx_fio_fopen(infile,"r");
1320 fgets2(buf,STRLEN-1,fp);
1321 if (sscanf(buf,"%d",&n) != 1)
1322 gmx_fatal(FARGS,"Can not read number of atoms from %s",infile);
1323 fgets2(buf,STRLEN-1,fp);
1324 strcpy(title,buf);
1325 for(i=0; (i<n); i++) {
1326 fgets2(buf,STRLEN-1,fp);
1327 if (sscanf(buf,"%s%lf%lf%lf",atomnm,&xx,&yy,&zz) != 4)
1328 gmx_fatal(FARGS,"Can not read coordinates from %s",infile);
1329 atoms->atomname[i] = put_symtab(tab,atomnm);
1330 x[i][XX] = xx*0.1;
1331 x[i][YY] = yy*0.1;
1332 x[i][ZZ] = zz*0.1;
1333 }
1334 gmx_fio_fclose(fp);
1335 }
1336
1337 void get_stx_coordnum(const char *infile,int *natoms)
1338 {
1339 FILE *in;
1340 int ftp,tpxver,tpxgen;
1341 t_trxframe fr;
1342
1343 ftp=fn2ftp(infile);
1344 range_check(ftp,0,efNR);
1345 switch (ftp) {
1346 case efGRO:
1347 get_coordnum(infile, natoms);
1348 break;
1349 case efG96:
1350 in=gmx_fio_fopen(infile,"r");
1351 fr.title = NULL;
1352 fr.natoms = -1;
1353 fr.atoms = NULL;
1354 fr.x = NULL;
1355 fr.v = NULL;
1356 fr.f = NULL;
1357 *natoms=read_g96_conf(in,infile,&fr);
1358 gmx_fio_fclose(in);
1359 break;
1360 case efXYZ:
1361 *natoms = get_xyz_coordnum(infile);
1362 break;
1363 case efPDB:
1364 case efBRK:
1365 case efENT:
1366 in=gmx_fio_fopen(infile,"r");
1367 get_pdb_coordnum(in, natoms);
1368 gmx_fio_fclose(in);
1369 break;
1370 case efESP:
1371 *natoms = get_espresso_coordnum(infile);
1372 break;
1373 case efTPA:
1374 case efTPB:
1375 case efTPR: {
1376 t_tpxheader tpx;
1377
1378 read_tpxheader(infile,&tpx,TRUE,&tpxver,&tpxgen);
1379 *natoms = tpx.natoms;
1380 break;
1381 }
1382 default:
1383 gmx_fatal(FARGS,"File type %s not supported in get_stx_coordnum",
1384 ftp2ext(ftp));
1385 }
1386 }
1387
1388 void read_stx_conf(const char *infile,char *title,t_atoms *atoms,
1389 rvec x[],rvec *v,int *ePBC,matrix box)
1390 {
1391 FILE *in;
1392 char buf[256];
1393 gmx_mtop_t *mtop;
1394 t_topology top;
1395 t_trxframe fr;
1396 int i,ftp,natoms;
1397 real d;
1398
1399 if (atoms->nr == 0)
1400 fprintf(stderr,"Warning: Number of atoms in %s is 0\n",infile);
1401 else if (atoms->atom == NULL)
1402 gmx_mem("Uninitialized array atom");
1403
1404 if (ePBC)
1405 *ePBC = -1;
1406
1407 ftp=fn2ftp(infile);
1408 switch (ftp) {
1409 case efGRO:
1410 read_whole_conf(infile, title, atoms, x, v, box);
1411 break;
1412 case efXYZ:
1413 read_xyz_conf(infile,title,atoms,x);
1414 break;
1415 case efG96:
1416 fr.title = title;
1417 fr.natoms = atoms->nr;
1418 fr.atoms = atoms;
1419 fr.x = x;
1420 fr.v = v;
1421 fr.f = NULL;
1422 in = gmx_fio_fopen(infile,"r");
1423 read_g96_conf(in, infile, &fr);
1424 gmx_fio_fclose(in);
1425 copy_mat(fr.box,box);
1426 break;
1427 case efPDB:
1428 case efBRK:
1429 case efENT:
1430 read_pdb_conf(infile, title, atoms, x, ePBC, box, TRUE, NULL);
1431 break;
1432 case efESP:
1433 read_espresso_conf(infile,atoms,x,v,box);
1434 break;
1435 case efTPR:
1436 case efTPB:
1437 case efTPA:
1438 snew(mtop,1);
1439 i = read_tpx(infile,NULL,box,&natoms,x,v,NULL,mtop);
1440 if (ePBC)
1441 *ePBC = i;
1442
1443 strcpy(title,*(mtop->name));
1444
1445 /* Free possibly allocated memory */
1446 done_atom(atoms);
1447
1448 *atoms = gmx_mtop_global_atoms(mtop);
1449 top = gmx_mtop_t_to_t_topology(mtop);
1450 tpx_make_chain_identifiers(atoms,&top.mols);
1451
1452 sfree(mtop);
1453 /* The strings in the symtab are still in use in the returned t_atoms
1454 * structure, so we should not free them. But there is no place to put the
1455 * symbols; the only choice is to leak the memory...
1456 * So we clear the symbol table before freeing the topology structure. */
1457 open_symtab(&top.symtab);
1458 done_top(&top);
1459
1460 break;
1461 default:
1462 gmx_incons("Not supported in read_stx_conf");
1463 }
1464 }
1465