3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
33 * GROningen Mixture of Alchemy and Childrens' Stories
55 #include "gmx_fatal.h"
60 #include "mtop_util.h"
65 static int read_g96_pos(char line[],t_symtab *symtab,
66 FILE *fp,const char *infile,
71 int nwanted,natoms,atnr,resnr=0,oldres,newres,shift;
72 char anm[STRLEN],resnm[STRLEN];
88 oldres = -666; /* Unlikely number for the first residue! */
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",
96 if ((nwanted != -1) && (natoms >= nwanted))
98 "Found more coordinates (%d) in %s than expected %d\n",
99 natoms,infile,nwanted);
102 (sscanf(line,"%5d%c%5s%c%5s%7d",&resnr,&c1,resnm,&c2,anm,&atnr)
108 strncpy(resnm,"???",sizeof(resnm)-1);
110 strncpy(anm,"???",sizeof(anm)-1);
112 atoms->atomname[natoms]=put_symtab(symtab,anm);
113 if (resnr != oldres) {
116 if (newres >= atoms->nr)
117 gmx_fatal(FARGS,"More residues than atoms in %s (natoms = %d)",
119 atoms->atom[natoms].resind = newres;
120 if (newres+1 > atoms->nres) {
121 atoms->nres = newres+1;
123 t_atoms_set_resinfo(atoms,natoms,symtab,resnm,resnr,' ',0,' ');
125 atoms->atom[natoms].resind = newres;
129 fr->x[natoms][0] = db1;
130 fr->x[natoms][1] = db2;
131 fr->x[natoms][2] = db3;
136 if ((nwanted != -1) && natoms != nwanted)
138 "Warning: found less coordinates (%d) in %s than expected %d\n",
139 natoms,infile,nwanted);
147 static int read_g96_vel(char line[],FILE *fp,const char *infile,
151 int nwanted,natoms=-1,shift;
154 nwanted = fr->natoms;
156 if (fr->v && fr->bV) {
157 if (strcmp(line,"VELOCITYRED") == 0)
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",
169 if ((nwanted != -1) && (natoms >= nwanted))
170 gmx_fatal(FARGS,"Found more velocities (%d) in %s than expected %d\n",
171 natoms,infile,nwanted);
173 fr->v[natoms][0] = db1;
174 fr->v[natoms][1] = db2;
175 fr->v[natoms][2] = db3;
180 if ((nwanted != -1) && (natoms != nwanted))
182 "Warning: found less velocities (%d) in %s than expected %d\n",
183 natoms,infile,nwanted);
189 int read_g96_conf(FILE *fp,const char *infile,t_trxframe *fr, char *line)
191 t_symtab *symtab=NULL;
192 gmx_bool bAtStart,bTime,bAtoms,bPos,bVel,bBox,bEnd,bFinished;
194 double db1,db2,db3,db4,db5,db6,db7,db8,db9;
196 bAtStart = (ftell(fp) == 0);
198 clear_trxframe(fr,FALSE);
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);
215 while (!bEnd && fgets2(line,STRLEN,fp))
216 bEnd = (strcmp(line,"END") == 0);
217 fgets2(line,STRLEN,fp);
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). */
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);
232 if (!fr->bTime && !fr->bX) {
236 bFinished = (fgets2(line,STRLEN,fp) == NULL);
237 while (!bFinished && (line[0] == '#'));
238 sscanf(line,"%15d%15lf",&(fr->step),&db1);
247 natoms = read_g96_pos(line,symtab,fp,infile,fr);
253 natoms = read_g96_vel(line,fp,infile,fr);
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);
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;
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;
281 } while (!bFinished && fgets2(line,STRLEN,fp));
290 void write_g96_conf(FILE *out,t_trxframe *fr,
291 int nindex,atom_id *index)
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);
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]);
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]);
327 fprintf(out,"END\n");
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]);
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]);
348 fprintf(out,"END\n");
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]);
360 fprintf(out,"END\n");
364 static int get_espresso_word(FILE *fp,char word[])
374 if (i == ' ' || i == '\n' || i == '\t') {
377 } else if (i == '{') {
381 } else if (i == '}') {
386 word[nc++] = (char)i;
389 } while (i != EOF && ret == 0);
393 /* printf("'%s'\n",word); */
398 static int check_open_parenthesis(FILE *fp,int r,
399 const char *infile,const char *keyword)
408 r = get_espresso_word(fp,word);
412 gmx_fatal(FARGS,"Expected '{' after '%s' in file '%s'",
419 static int check_close_parenthesis(FILE *fp,int r,
420 const char *infile,const char *keyword)
429 r = get_espresso_word(fp,word);
433 gmx_fatal(FARGS,"Expected '}' after section '%s' in file '%s'",
440 enum { espID, espPOS, espTYPE, espQ, espV, espF, espMOLECULE, espNR };
441 const char *esp_prop[espNR] = { "id", "pos", "type", "q", "v", "f",
444 static void read_espresso_conf(const char *infile,
445 t_atoms *atoms,rvec x[],rvec *v,matrix box)
447 t_symtab *symtab=NULL;
449 char word[STRLEN],buf[STRLEN];
450 int natoms,level,npar,r,nprop,p,i,m,molnr;
453 gmx_bool bFoundParticles,bFoundProp,bFoundVariable,bMol;
462 fp = gmx_fio_fopen(infile,"r");
464 bFoundParticles = FALSE;
465 bFoundVariable = FALSE;
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");
473 while (level == 2 && (r=get_espresso_word(fp,word))) {
475 for(p=0; p<espNR; p++) {
476 if (strcmp(word,esp_prop[p]) == 0) {
479 /* printf(" prop[%d] = %s\n",nprop-1,esp_prop[prop[nprop-1]]); */
482 if (!bFoundProp && word[0] != '}') {
483 gmx_fatal(FARGS,"Can not read Espresso files with particle property '%s'",word);
485 if (bFoundProp && p == espMOLECULE)
492 while (level > 0 && (r=get_espresso_word(fp,word))) {
499 for(p=0; p<nprop; p++) {
502 r = get_espresso_word(fp,word);
507 r = get_espresso_word(fp,word);
508 sscanf(word,"%lf",&d);
513 r = get_espresso_word(fp,word);
514 atoms->atom[i].type = strtol(word, NULL, 10);
517 r = get_espresso_word(fp,word);
518 sscanf(word,"%lf",&d);
519 atoms->atom[i].q = d;
523 r = get_espresso_word(fp,word);
524 sscanf(word,"%lf",&d);
530 r = get_espresso_word(fp,word);
535 r = get_espresso_word(fp,word);
536 molnr = strtol(word, NULL, 10);
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? */
546 atoms->atom[i].resind = atoms->atom[i-1].resind;
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);
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");
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);
567 'A'+atoms->atom[i].type/26,'A'+atoms->atom[i].type%26);
569 t_atoms_set_resinfo(atoms,i,symtab,buf,i,' ',0,' ');
577 atoms->nres = atoms->nr;
579 if (i != atoms->nr) {
580 gmx_fatal(FARGS,"Internal inconsistency in Espresso routines, read %d atoms, expected %d atoms",i,atoms->nr);
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) {
588 r = get_espresso_word(fp,word);
589 sscanf(word,"%lf",&d);
592 level += check_close_parenthesis(fp,r,infile,"box_l");
602 if (!bFoundParticles) {
603 fprintf(stderr,"Did not find a particles section in Espresso file '%s'\n",
610 static int get_espresso_coordnum(const char *infile)
615 gmx_bool bFoundParticles;
619 fp = gmx_fio_fopen(infile,"r");
621 bFoundParticles = FALSE;
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))) {
642 if (!bFoundParticles) {
643 fprintf(stderr,"Did not find a particles section in Espresso file '%s'\n",
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)
658 fprintf(out,"# %s\n",title);
659 if (TRICLINIC(box)) {
660 gmx_warning("The Espresso format does not support triclinic unit-cells");
662 fprintf(out,"{variable {box_l %f %f %f}}\n",box[0][0],box[1][1],box[2][2]);
664 fprintf(out,"{particles {id pos type q%s}\n",v ? " v" : "");
665 for(i=0; i<nx; 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);
674 fprintf(out," %f %f %f",v[j][XX],v[j][YY],v[j][ZZ]);
680 static void get_coordnum_fp (FILE *in, char *title, int *natoms)
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");
691 static void get_coordnum (const char *infile,int *natoms)
696 in=gmx_fio_fopen(infile,"r");
697 get_coordnum_fp(in,title,natoms);
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)
706 char line[STRLEN+1],*ptr;
708 double x1,y1,z1,x2,y2,z2;
710 int natoms,i,m,resnr,newres,oldres,ddist,c;
711 gmx_bool bFirst,bVel;
715 oldres = NOTSET; /* Unlikely number for the first residue! */
718 /* Read the title and number of atoms */
719 get_coordnum_fp(in,title,&natoms);
721 if (natoms > atoms->nr)
722 gmx_fatal(FARGS,"gro file contains more atoms (%d) than expected (%d)",
724 else if (natoms < atoms->nr)
725 fprintf(stderr,"Warning: gro file contains less atoms (%d) than expected"
726 " (%d)\n",natoms,atoms->nr);
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);
737 if (strlen(line) < 39)
738 gmx_fatal(FARGS,"Invalid line in %s for atom %d:\n%s",infile,i+1,line);
740 /* determine read precision from distance between periods
746 gmx_fatal(FARGS,"A coordinate in file %s does not contain a '.'",infile);
747 p2=strchr(&p1[1],'.');
749 gmx_fatal(FARGS,"A coordinate in file %s does not contain a '.'",infile);
753 p3=strchr(&p2[1],'.');
755 gmx_fatal(FARGS,"A coordinate in file %s does not contain a '.'",infile);
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);
764 sscanf(name,"%d",&resnr);
765 memcpy(name,line+5,5);
767 if (resnr != oldres) {
770 if (newres >= natoms)
771 gmx_fatal(FARGS,"More residues than atoms in %s (natoms = %d)",
773 atoms->atom[i].resind = newres;
774 t_atoms_set_resinfo(atoms,i,symtab,name,resnr,' ',0,' ');
776 atoms->atom[i].resind = newres;
780 memcpy(name,line+10,5);
781 atoms->atomname[i]=put_symtab(symtab,name);
783 /* eventueel controle atomnumber met i+1 */
785 /* coordinates (start after residue shit) */
787 /* Read fixed format */
788 for(m=0; m<DIM; m++) {
789 for(c=0; (c<ddist && ptr[0]); c++) {
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);
801 /* velocities (start after residues and coordinates) */
803 /* Read fixed format */
804 for(m=0; m<DIM; m++) {
805 for(c=0; (c<ddist && ptr[0]); c++) {
810 if (sscanf (buf,"%lf",&x1) != 1) {
819 atoms->nres = newres + 1;
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);
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]);
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]);
841 /* We found the first three values, the diagonal elements */
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;
859 static void read_whole_conf(const char *infile,char *title,
860 t_atoms *atoms, rvec x[],rvec *v, matrix box)
867 in=gmx_fio_fopen(infile,"r");
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);
878 gmx_bool gro_next_x_or_v(FILE *status,t_trxframe *fr)
882 char title[STRLEN],*p;
889 open_symtab(&symtab);
891 snew(atoms.atom,fr->natoms);
892 atoms.nres=fr->natoms;
893 snew(atoms.resinfo,fr->natoms);
894 snew(atoms.atomname,fr->natoms);
896 fr->bV = get_w_conf(status,title,title,&symtab,&atoms,&ndec,fr->x,fr->v,fr->box);
899 /* prec = 10^ndec: */
900 for(i=0; i<ndec; i++)
908 sfree(atoms.resinfo);
909 sfree(atoms.atomname);
910 done_symtab(&symtab);
912 if ((p=strstr(title,"t=")) != NULL) {
914 if (sscanf(p,"%lf",&tt)==1) {
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);
929 int gro_first_x_or_v(FILE *status,t_trxframe *fr)
934 fprintf(stderr,"Reading frames from gro file");
935 get_coordnum_fp(status, title, &fr->natoms);
937 fprintf(stderr," '%s', %d atoms.\n",title, fr->natoms);
941 gmx_file("No coordinates in gro file");
943 snew(fr->x,fr->natoms);
944 snew(fr->v,fr->natoms);
945 gro_next_x_or_v(status, fr);
950 static void make_hconf_format(int pr,gmx_bool bVel,char format[])
954 /* build format string for printing,
955 something like "%8.3f" for x and "%8.4f" for v */
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);
966 sprintf(format,"%%%d.%df%%%d.%df%%%d.%df\n",l,pr,l,pr,l,pr);
970 static void write_hconf_box(FILE *out,int pr,matrix box)
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);
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]);
989 sprintf(format,"%%%d.%df%%%d.%df%%%d.%df\n",l,pr,l,pr,l,pr);
991 box[XX][XX],box[YY][YY],box[ZZ][ZZ]);
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)
999 char resnm[6],nm[6],format[100];
1000 int ai,i,resind,resnr;
1003 fprintf (out,"%s\n",(title && title[0])?title:format);
1004 fprintf (out,"%5d\n",nx);
1006 make_hconf_format(pr,v!=NULL,format);
1008 for (i=0; (i<nx); i++) {
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;
1017 strncpy(resnm," ??? ",sizeof(resnm)-1);
1022 strncpy(nm,*atoms->atomname[ai],sizeof(nm)-1);
1024 strncpy(nm," ??? ",sizeof(nm)-1);
1026 fprintf(out,"%5d%-5.5s%5.5s%5d",resnr%100000,resnm,nm,(ai+1)%100000);
1027 /* next fprintf uses built format string */
1030 x[ai][XX], x[ai][YY], x[ai][ZZ], v[ai][XX],v[ai][YY],v[ai][ZZ]);
1033 x[ai][XX], x[ai][YY], x[ai][ZZ]);
1036 write_hconf_box(out,pr,box);
1041 static void write_hconf_mtop(FILE *out,const char *title,gmx_mtop_t *mtop,
1043 rvec *x,rvec *v,matrix box)
1047 gmx_mtop_atomloop_all_t aloop;
1049 char *atomname,*resname;
1052 fprintf (out,"%s\n",(title && title[0])?title:format);
1053 fprintf (out,"%5d\n",mtop->natoms);
1055 make_hconf_format(pr,v!=NULL,format);
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);
1061 fprintf(out,"%5d%-5.5s%5.5s%5d",
1062 resnr%100000,resname,atomname,(i+1)%100000);
1063 /* next fprintf uses built format string */
1066 x[i][XX], x[i][YY], x[i][ZZ], v[i][XX],v[i][YY],v[i][ZZ]);
1069 x[i][XX], x[i][YY], x[i][ZZ]);
1072 write_hconf_box(out,pr,box);
1077 void write_hconf_p(FILE *out,const char *title,t_atoms *atoms, int pr,
1078 rvec *x,rvec *v,matrix box)
1084 for(i=0; (i<atoms->nr); i++)
1086 write_hconf_indexed_p(out,title,atoms,atoms->nr,aa,pr,x,v,box);
1090 void write_conf_p(const char *outfile, const char *title,
1091 t_atoms *atoms, int pr,
1092 rvec *x, rvec *v,matrix box)
1096 out=gmx_fio_fopen(outfile,"w");
1097 write_hconf_p(out,title,atoms,pr,x,v,box);
1099 gmx_fio_fclose (out);
1102 static void write_conf(const char *outfile, const char *title, t_atoms *atoms,
1103 rvec *x, rvec *v,matrix box)
1105 write_conf_p(outfile, title, atoms, 3, x, v, box);
1108 void write_sto_conf_indexed(const char *outfile,const char *title,
1110 rvec x[],rvec *v,int ePBC,matrix box,
1111 atom_id nindex,atom_id index[])
1117 ftp=fn2ftp(outfile);
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);
1125 clear_trxframe(&fr,TRUE);
1128 fr.natoms = atoms->nr;
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);
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);
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);
1159 gmx_fatal(FARGS,"Sorry, can not write a topology to %s",outfile);
1162 gmx_incons("Not supported in write_sto_conf_indexed");
1166 static void write_xyz_conf(const char *outfile,const char *title,
1167 t_atoms *atoms,rvec *x)
1173 gmx_atomprop_t aps=gmx_atomprop_init();
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);
1185 if ((ptr = gmx_atomprop_element(aps,anr)) == NULL)
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]);
1191 gmx_atomprop_destroy(aps);
1194 void write_sto_conf(const char *outfile,const char *title,t_atoms *atoms,
1195 rvec x[],rvec *v,int ePBC,matrix box)
1201 ftp=fn2ftp(outfile);
1204 write_conf(outfile, title, atoms, x, v, box);
1207 clear_trxframe(&fr,TRUE);
1210 fr.natoms = atoms->nr;
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);
1226 write_xyz_conf(outfile,(strlen(title) > 0) ? title : outfile,atoms,x);
1231 out=gmx_fio_fopen(outfile,"w");
1232 write_pdbfile(out, title, atoms, x, ePBC, box, ' ', -1,NULL,TRUE);
1233 gmx_fio_fclose(out);
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);
1243 gmx_fatal(FARGS,"Sorry, can not write a topology to %s",outfile);
1246 gmx_incons("Not supported in write_sto_conf");
1250 void write_sto_conf_mtop(const char *outfile,const char *title,
1252 rvec x[],rvec *v,int ePBC,matrix box)
1258 ftp=fn2ftp(outfile);
1261 out = gmx_fio_fopen(outfile,"w");
1262 write_hconf_mtop(out,title,mtop,3,x,v,box);
1263 gmx_fio_fclose(out);
1266 /* This is a brute force approach which requires a lot of memory.
1267 * We should implement mtop versions of all writing routines.
1269 atoms = gmx_mtop_global_atoms(mtop);
1271 write_sto_conf(outfile,title,&atoms,x,v,ePBC,box);
1278 static int get_xyz_coordnum(const char *infile)
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);
1291 static void read_xyz_conf(const char *infile,char *title,
1292 t_atoms *atoms,rvec *x)
1298 char atomnm[32],buf[STRLEN];
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);
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);
1319 void get_stx_coordnum(const char *infile,int *natoms)
1322 int ftp,tpxver,tpxgen;
1324 char g96_line[STRLEN+1];
1327 range_check(ftp,0,efNR);
1330 get_coordnum(infile, natoms);
1333 in=gmx_fio_fopen(infile,"r");
1340 *natoms=read_g96_conf(in,infile,&fr,g96_line);
1344 *natoms = get_xyz_coordnum(infile);
1349 in=gmx_fio_fopen(infile,"r");
1350 get_pdb_coordnum(in, natoms);
1354 *natoms = get_espresso_coordnum(infile);
1361 read_tpxheader(infile,&tpx,TRUE,&tpxver,&tpxgen);
1362 *natoms = tpx.natoms;
1366 gmx_fatal(FARGS,"File type %s not supported in get_stx_coordnum",
1371 void read_stx_conf(const char *infile,char *title,t_atoms *atoms,
1372 rvec x[],rvec *v,int *ePBC,matrix box)
1381 char g96_line[STRLEN+1];
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");
1394 read_whole_conf(infile, title, atoms, x, v, box);
1397 read_xyz_conf(infile,title,atoms,x);
1401 fr.natoms = atoms->nr;
1406 in = gmx_fio_fopen(infile,"r");
1407 read_g96_conf(in, infile, &fr, g96_line);
1409 copy_mat(fr.box,box);
1414 read_pdb_conf(infile, title, atoms, x, ePBC, box, TRUE, NULL);
1417 read_espresso_conf(infile,atoms,x,v,box);
1423 i = read_tpx(infile,NULL,box,&natoms,x,v,NULL,mtop);
1427 strcpy(title,*(mtop->name));
1429 /* Free possibly allocated memory */
1432 *atoms = gmx_mtop_global_atoms(mtop);
1433 top = gmx_mtop_t_to_t_topology(mtop);
1434 tpx_make_chain_identifiers(atoms,&top.mols);
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);
1446 gmx_incons("Not supported in read_stx_conf");