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,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);
101 if (atoms && fr->bAtoms &&
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)
191 t_symtab *symtab=NULL;
193 bool bAtStart,bTime,bAtoms,bPos,bVel,bBox,bEnd,bFinished;
195 double db1,db2,db3,db4,db5,db6,db7,db8,db9;
197 bAtStart = (ftell(fp) == 0);
199 clear_trxframe(fr,FALSE);
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);
216 while (!bEnd && fgets2(line,STRLEN,fp))
217 bEnd = (strcmp(line,"END") == 0);
218 fgets2(line,STRLEN,fp);
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). */
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);
233 if (!fr->bTime && !fr->bX) {
237 bFinished = (fgets2(line,STRLEN,fp) == NULL);
238 while (!bFinished && (line[0] == '#'));
239 sscanf(line,"%15d%15lf",&(fr->step),&db1);
248 natoms = read_g96_pos(line,symtab,fp,infile,fr);
254 natoms = read_g96_vel(line,fp,infile,fr);
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);
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;
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;
282 } while (!bFinished && fgets2(line,STRLEN,fp));
291 void write_g96_conf(FILE *out,t_trxframe *fr,
292 int nindex,atom_id *index)
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);
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]);
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]);
328 fprintf(out,"END\n");
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]);
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]);
349 fprintf(out,"END\n");
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]);
361 fprintf(out,"END\n");
365 static int get_espresso_word(FILE *fp,char word[])
375 if (i == ' ' || i == '\n' || i == '\t') {
378 } else if (i == '{') {
382 } else if (i == '}') {
387 word[nc++] = (char)i;
390 } while (i != EOF && ret == 0);
394 /* printf("'%s'\n",word); */
399 static int check_open_parenthesis(FILE *fp,int r,
400 const char *infile,const char *keyword)
409 r = get_espresso_word(fp,word);
413 gmx_fatal(FARGS,"Expected '{' after '%s' in file '%s'",
420 static int check_close_parenthesis(FILE *fp,int r,
421 const char *infile,const char *keyword)
430 r = get_espresso_word(fp,word);
434 gmx_fatal(FARGS,"Expected '}' after section '%s' in file '%s'",
441 enum { espID, espPOS, espTYPE, espQ, espV, espF, espMOLECULE, espNR };
442 const char *esp_prop[espNR] = { "id", "pos", "type", "q", "v", "f",
445 static void read_espresso_conf(const char *infile,
446 t_atoms *atoms,rvec x[],rvec *v,matrix box)
448 t_symtab *symtab=NULL;
450 char word[STRLEN],buf[STRLEN];
451 int natoms,level,npar,r,nprop,p,i,m,molnr;
454 bool bFoundParticles,bFoundProp,bFoundVariable,bMol;
463 fp = gmx_fio_fopen(infile,"r");
465 bFoundParticles = FALSE;
466 bFoundVariable = FALSE;
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");
474 while (level == 2 && (r=get_espresso_word(fp,word))) {
476 for(p=0; p<espNR; p++) {
477 if (strcmp(word,esp_prop[p]) == 0) {
480 /* printf(" prop[%d] = %s\n",nprop-1,esp_prop[prop[nprop-1]]); */
483 if (!bFoundProp && word[0] != '}') {
484 gmx_fatal(FARGS,"Can not read Espresso files with particle property '%s'",word);
486 if (bFoundProp && p == espMOLECULE)
493 while (level > 0 && (r=get_espresso_word(fp,word))) {
500 for(p=0; p<nprop; p++) {
503 r = get_espresso_word(fp,word);
508 r = get_espresso_word(fp,word);
509 sscanf(word,"%lf",&d);
514 r = get_espresso_word(fp,word);
515 atoms->atom[i].type = strtol(word, NULL, 10);
518 r = get_espresso_word(fp,word);
519 sscanf(word,"%lf",&d);
520 atoms->atom[i].q = d;
524 r = get_espresso_word(fp,word);
525 sscanf(word,"%lf",&d);
531 r = get_espresso_word(fp,word);
536 r = get_espresso_word(fp,word);
537 molnr = strtol(word, NULL, 10);
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? */
547 atoms->atom[i].resind = atoms->atom[i-1].resind;
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);
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");
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);
568 'A'+atoms->atom[i].type/26,'A'+atoms->atom[i].type%26);
570 t_atoms_set_resinfo(atoms,i,symtab,buf,i,' ',0,' ');
578 atoms->nres = atoms->nr;
580 if (i != atoms->nr) {
581 gmx_fatal(FARGS,"Internal inconsistency in Espresso routines, read %d atoms, expected %d atoms",i,atoms->nr);
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) {
589 r = get_espresso_word(fp,word);
590 sscanf(word,"%lf",&d);
593 level += check_close_parenthesis(fp,r,infile,"box_l");
603 if (!bFoundParticles) {
604 fprintf(stderr,"Did not find a particles section in Espresso file '%s'\n",
611 static int get_espresso_coordnum(const char *infile)
616 bool bFoundParticles;
620 fp = gmx_fio_fopen(infile,"r");
622 bFoundParticles = FALSE;
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))) {
643 if (!bFoundParticles) {
644 fprintf(stderr,"Did not find a particles section in Espresso file '%s'\n",
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)
659 fprintf(out,"# %s\n",title);
660 if (TRICLINIC(box)) {
661 gmx_warning("The Espresso format does not support triclinic unit-cells");
663 fprintf(out,"{variable {box_l %f %f %f}}\n",box[0][0],box[1][1],box[2][2]);
665 fprintf(out,"{particles {id pos type q%s}\n",v ? " v" : "");
666 for(i=0; i<nx; 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);
675 fprintf(out," %f %f %f",v[j][XX],v[j][YY],v[j][ZZ]);
681 static void get_coordnum_fp (FILE *in, char *title, int *natoms)
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");
692 static void get_coordnum (const char *infile,int *natoms)
697 in=gmx_fio_fopen(infile,"r");
698 get_coordnum_fp(in,title,natoms);
702 static bool get_w_conf(FILE *in,const char *infile,char *title,
703 t_atoms *atoms, int *ndec, rvec x[],rvec *v, matrix box)
705 t_symtab *symtab=NULL;
707 char line[STRLEN+1],*ptr;
709 double x1,y1,z1,x2,y2,z2;
711 int natoms,i,m,resnr,newres,oldres,ddist,c;
716 oldres = NOTSET; /* Unlikely number for the first residue! */
724 /* Read the title and number of atoms */
725 get_coordnum_fp(in,title,&natoms);
727 if (natoms > atoms->nr)
728 gmx_fatal(FARGS,"gro file contains more atoms (%d) than expected (%d)",
730 else if (natoms < atoms->nr)
731 fprintf(stderr,"Warning: gro file contains less atoms (%d) than expected"
732 " (%d)\n",natoms,atoms->nr);
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);
743 if (strlen(line) < 39)
744 gmx_fatal(FARGS,"Invalid line in %s for atom %d:\n%s",infile,i+1,line);
746 /* determine read precision from distance between periods
752 gmx_fatal(FARGS,"A coordinate in file %s does not contain a '.'",infile);
753 p2=strchr(&p1[1],'.');
755 gmx_fatal(FARGS,"A coordinate in file %s does not contain a '.'",infile);
759 p3=strchr(&p2[1],'.');
761 gmx_fatal(FARGS,"A coordinate in file %s does not contain a '.'",infile);
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);
770 sscanf(name,"%d",&resnr);
771 memcpy(name,line+5,5);
773 if (resnr != oldres) {
776 if (newres >= natoms)
777 gmx_fatal(FARGS,"More residues than atoms in %s (natoms = %d)",
779 atoms->atom[i].resind = newres;
780 t_atoms_set_resinfo(atoms,i,symtab,name,resnr,' ',0,' ');
782 atoms->atom[i].resind = newres;
786 memcpy(name,line+10,5);
787 atoms->atomname[i]=put_symtab(symtab,name);
789 /* eventueel controle atomnumber met i+1 */
791 /* coordinates (start after residue shit) */
793 /* Read fixed format */
794 for(m=0; m<DIM; m++) {
795 for(c=0; (c<ddist && ptr[0]); c++) {
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);
807 /* velocities (start after residues and coordinates) */
809 /* Read fixed format */
810 for(m=0; m<DIM; m++) {
811 for(c=0; (c<ddist && ptr[0]); c++) {
816 if (sscanf (buf,"%lf",&x1) != 1) {
825 atoms->nres = newres + 1;
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);
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]);
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]);
847 /* We found the first three values, the diagonal elements */
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;
861 close_symtab(symtab);
866 static void read_whole_conf(const char *infile,char *title,
867 t_atoms *atoms, rvec x[],rvec *v, matrix box)
873 in=gmx_fio_fopen(infile,"r");
875 get_w_conf(in, infile, title, atoms, &ndec, x, v, box);
880 static void get_conf(FILE *in,char *title,int *natoms,
881 rvec x[],rvec *v,matrix box)
887 snew(atoms.atom,*natoms);
889 snew(atoms.resinfo,*natoms);
890 snew(atoms.atomname,*natoms);
892 get_w_conf(in,title,title,&atoms,&ndec,x,v,box);
895 sfree(atoms.resinfo);
896 sfree(atoms.atomname);
899 bool gro_next_x_or_v(FILE *status,t_trxframe *fr)
902 char title[STRLEN],*p;
910 snew(atoms.atom,fr->natoms);
911 atoms.nres=fr->natoms;
912 snew(atoms.resinfo,fr->natoms);
913 snew(atoms.atomname,fr->natoms);
915 fr->bV = get_w_conf(status,title,title,&atoms,&ndec,fr->x,fr->v,fr->box);
918 /* prec = 10^ndec: */
919 for(i=0; i<ndec; i++)
927 sfree(atoms.resinfo);
928 sfree(atoms.atomname);
930 if ((p=strstr(title,"t=")) != NULL) {
932 if (sscanf(p,"%lf",&tt)==1) {
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);
947 int gro_first_x_or_v(FILE *status,t_trxframe *fr)
952 fprintf(stderr,"Reading frames from gro file");
953 get_coordnum_fp(status, title, &fr->natoms);
955 fprintf(stderr," '%s', %d atoms.\n",title, fr->natoms);
959 gmx_file("No coordinates in gro file");
961 snew(fr->x,fr->natoms);
962 snew(fr->v,fr->natoms);
963 gro_next_x_or_v(status, fr);
968 static void make_hconf_format(int pr,bool bVel,char format[])
972 /* build format string for printing,
973 something like "%8.3f" for x and "%8.4f" for v */
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);
984 sprintf(format,"%%%d.%df%%%d.%df%%%d.%df\n",l,pr,l,pr,l,pr);
988 static void write_hconf_box(FILE *out,int pr,matrix box)
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);
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]);
1007 sprintf(format,"%%%d.%df%%%d.%df%%%d.%df\n",l,pr,l,pr,l,pr);
1009 box[XX][XX],box[YY][YY],box[ZZ][ZZ]);
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)
1017 char resnm[6],nm[6],format[100];
1018 int ai,i,resind,resnr;
1021 fprintf (out,"%s\n",(title && title[0])?title:format);
1022 fprintf (out,"%5d\n",nx);
1024 make_hconf_format(pr,v!=NULL,format);
1026 for (i=0; (i<nx); i++) {
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;
1035 strncpy(resnm," ??? ",sizeof(resnm)-1);
1040 strncpy(nm,*atoms->atomname[ai],sizeof(nm)-1);
1042 strncpy(nm," ??? ",sizeof(nm)-1);
1044 fprintf(out,"%5d%-5.5s%5.5s%5d",resnr%100000,resnm,nm,(ai+1)%100000);
1045 /* next fprintf uses built format string */
1048 x[ai][XX], x[ai][YY], x[ai][ZZ], v[ai][XX],v[ai][YY],v[ai][ZZ]);
1051 x[ai][XX], x[ai][YY], x[ai][ZZ]);
1054 write_hconf_box(out,pr,box);
1059 static void write_hconf_mtop(FILE *out,const char *title,gmx_mtop_t *mtop,
1061 rvec *x,rvec *v,matrix box)
1065 gmx_mtop_atomloop_all_t aloop;
1067 char *atomname,*resname;
1070 fprintf (out,"%s\n",(title && title[0])?title:format);
1071 fprintf (out,"%5d\n",mtop->natoms);
1073 make_hconf_format(pr,v!=NULL,format);
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);
1079 fprintf(out,"%5d%-5.5s%5.5s%5d",
1080 resnr%100000,resname,atomname,(i+1)%100000);
1081 /* next fprintf uses built format string */
1084 x[i][XX], x[i][YY], x[i][ZZ], v[i][XX],v[i][YY],v[i][ZZ]);
1087 x[i][XX], x[i][YY], x[i][ZZ]);
1090 write_hconf_box(out,pr,box);
1095 void write_hconf_p(FILE *out,const char *title,t_atoms *atoms, int pr,
1096 rvec *x,rvec *v,matrix box)
1102 for(i=0; (i<atoms->nr); i++)
1104 write_hconf_indexed_p(out,title,atoms,atoms->nr,aa,pr,x,v,box);
1108 void write_conf_p(const char *outfile, const char *title,
1109 t_atoms *atoms, int pr,
1110 rvec *x, rvec *v,matrix box)
1114 out=gmx_fio_fopen(outfile,"w");
1115 write_hconf_p(out,title,atoms,pr,x,v,box);
1117 gmx_fio_fclose (out);
1120 static void write_conf(const char *outfile, const char *title, t_atoms *atoms,
1121 rvec *x, rvec *v,matrix box)
1123 write_conf_p(outfile, title, atoms, 3, x, v, box);
1126 void write_sto_conf_indexed(const char *outfile,const char *title,
1128 rvec x[],rvec *v,int ePBC,matrix box,
1129 atom_id nindex,atom_id index[])
1135 ftp=fn2ftp(outfile);
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);
1143 clear_trxframe(&fr,TRUE);
1146 fr.natoms = atoms->nr;
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);
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);
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);
1177 gmx_fatal(FARGS,"Sorry, can not write a topology to %s",outfile);
1180 gmx_incons("Not supported in write_sto_conf_indexed");
1184 static void write_xyz_conf(const char *outfile,const char *title,
1185 t_atoms *atoms,rvec *x)
1191 gmx_atomprop_t aps=gmx_atomprop_init();
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);
1203 if ((ptr = gmx_atomprop_element(aps,anr)) == NULL)
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]);
1209 gmx_atomprop_destroy(aps);
1212 void write_sto_conf(const char *outfile,const char *title,t_atoms *atoms,
1213 rvec x[],rvec *v,int ePBC,matrix box)
1219 ftp=fn2ftp(outfile);
1222 write_conf(outfile, title, atoms, x, v, box);
1225 clear_trxframe(&fr,TRUE);
1228 fr.natoms = atoms->nr;
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);
1244 write_xyz_conf(outfile,(strlen(title) > 0) ? title : outfile,atoms,x);
1249 out=gmx_fio_fopen(outfile,"w");
1250 write_pdbfile(out, title, atoms, x, ePBC, box, ' ', -1,NULL,TRUE);
1251 gmx_fio_fclose(out);
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);
1261 gmx_fatal(FARGS,"Sorry, can not write a topology to %s",outfile);
1264 gmx_incons("Not supported in write_sto_conf");
1268 void write_sto_conf_mtop(const char *outfile,const char *title,
1270 rvec x[],rvec *v,int ePBC,matrix box)
1276 ftp=fn2ftp(outfile);
1279 out = gmx_fio_fopen(outfile,"w");
1280 write_hconf_mtop(out,title,mtop,3,x,v,box);
1281 gmx_fio_fclose(out);
1284 /* This is a brute force approach which requires a lot of memory.
1285 * We should implement mtop versions of all writing routines.
1287 atoms = gmx_mtop_global_atoms(mtop);
1289 write_sto_conf(outfile,title,&atoms,x,v,ePBC,box);
1296 static int get_xyz_coordnum(const char *infile)
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);
1309 static void read_xyz_conf(const char *infile,char *title,
1310 t_atoms *atoms,rvec *x)
1316 char atomnm[32],buf[STRLEN];
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);
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);
1337 void get_stx_coordnum(const char *infile,int *natoms)
1340 int ftp,tpxver,tpxgen;
1344 range_check(ftp,0,efNR);
1347 get_coordnum(infile, natoms);
1350 in=gmx_fio_fopen(infile,"r");
1357 *natoms=read_g96_conf(in,infile,&fr);
1361 *natoms = get_xyz_coordnum(infile);
1366 in=gmx_fio_fopen(infile,"r");
1367 get_pdb_coordnum(in, natoms);
1371 *natoms = get_espresso_coordnum(infile);
1378 read_tpxheader(infile,&tpx,TRUE,&tpxver,&tpxgen);
1379 *natoms = tpx.natoms;
1383 gmx_fatal(FARGS,"File type %s not supported in get_stx_coordnum",
1388 void read_stx_conf(const char *infile,char *title,t_atoms *atoms,
1389 rvec x[],rvec *v,int *ePBC,matrix box)
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");
1410 read_whole_conf(infile, title, atoms, x, v, box);
1413 read_xyz_conf(infile,title,atoms,x);
1417 fr.natoms = atoms->nr;
1422 in = gmx_fio_fopen(infile,"r");
1423 read_g96_conf(in, infile, &fr);
1425 copy_mat(fr.box,box);
1430 read_pdb_conf(infile, title, atoms, x, ePBC, box, TRUE, NULL);
1433 read_espresso_conf(infile,atoms,x,v,box);
1439 i = read_tpx(infile,NULL,box,&natoms,x,v,NULL,mtop);
1443 strcpy(title,*(mtop->name));
1445 /* Free possibly allocated memory */
1448 *atoms = gmx_mtop_global_atoms(mtop);
1449 top = gmx_mtop_t_to_t_topology(mtop);
1450 tpx_make_chain_identifiers(atoms,&top.mols);
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);
1462 gmx_incons("Not supported in read_stx_conf");