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