2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
43 #include "gromacs/utility/smalloc.h"
46 #include "gromacs/utility/cstringutil.h"
48 #include "gromacs/math/vec.h"
50 #include "gromacs/utility/futil.h"
56 #include "gromacs/utility/fatalerror.h"
59 #include "mtop_util.h"
64 static int read_g96_pos(char line
[], t_symtab
*symtab
,
65 FILE *fp
, const char *infile
,
70 int nwanted
, natoms
, atnr
, resnr
= 0, oldres
, newres
, shift
;
71 char anm
[STRLEN
], resnm
[STRLEN
];
92 oldres
= -666; /* Unlikely number for the first residue! */
94 while (!bEnd
&& fgets2(line
, STRLEN
, fp
))
96 bEnd
= (strncmp(line
, "END", 3) == 0);
97 if (!bEnd
&& (line
[0] != '#'))
99 if (sscanf(line
+shift
, "%15lf%15lf%15lf", &db1
, &db2
, &db3
) != 3)
101 gmx_fatal(FARGS
, "Did not find 3 coordinates for atom %d in %s\n",
104 if ((nwanted
!= -1) && (natoms
>= nwanted
))
107 "Found more coordinates (%d) in %s than expected %d\n",
108 natoms
, infile
, nwanted
);
113 (sscanf(line
, "%5d%c%5s%c%5s%7d", &resnr
, &c1
, resnm
, &c2
, anm
, &atnr
)
123 strncpy(resnm
, "???", sizeof(resnm
)-1);
125 strncpy(anm
, "???", sizeof(anm
)-1);
127 atoms
->atomname
[natoms
] = put_symtab(symtab
, anm
);
132 if (newres
>= atoms
->nr
)
134 gmx_fatal(FARGS
, "More residues than atoms in %s (natoms = %d)",
137 atoms
->atom
[natoms
].resind
= newres
;
138 if (newres
+1 > atoms
->nres
)
140 atoms
->nres
= newres
+1;
142 t_atoms_set_resinfo(atoms
, natoms
, symtab
, resnm
, resnr
, ' ', 0, ' ');
146 atoms
->atom
[natoms
].resind
= newres
;
151 fr
->x
[natoms
][0] = db1
;
152 fr
->x
[natoms
][1] = db2
;
153 fr
->x
[natoms
][2] = db3
;
158 if ((nwanted
!= -1) && natoms
!= nwanted
)
161 "Warning: found less coordinates (%d) in %s than expected %d\n",
162 natoms
, infile
, nwanted
);
171 static int read_g96_vel(char line
[], FILE *fp
, const char *infile
,
175 int nwanted
, natoms
= -1, shift
;
176 double db1
, db2
, db3
;
178 nwanted
= fr
->natoms
;
182 if (strcmp(line
, "VELOCITYRED") == 0)
192 while (!bEnd
&& fgets2(line
, STRLEN
, fp
))
194 bEnd
= (strncmp(line
, "END", 3) == 0);
195 if (!bEnd
&& (line
[0] != '#'))
197 if (sscanf(line
+shift
, "%15lf%15lf%15lf", &db1
, &db2
, &db3
) != 3)
199 gmx_fatal(FARGS
, "Did not find 3 velocities for atom %d in %s",
202 if ((nwanted
!= -1) && (natoms
>= nwanted
))
204 gmx_fatal(FARGS
, "Found more velocities (%d) in %s than expected %d\n",
205 natoms
, infile
, nwanted
);
209 fr
->v
[natoms
][0] = db1
;
210 fr
->v
[natoms
][1] = db2
;
211 fr
->v
[natoms
][2] = db3
;
216 if ((nwanted
!= -1) && (natoms
!= nwanted
))
219 "Warning: found less velocities (%d) in %s than expected %d\n",
220 natoms
, infile
, nwanted
);
227 int read_g96_conf(FILE *fp
, const char *infile
, t_trxframe
*fr
, char *line
)
229 t_symtab
*symtab
= NULL
;
230 gmx_bool bAtStart
, bTime
, bAtoms
, bPos
, bVel
, bBox
, bEnd
, bFinished
;
232 double db1
, db2
, db3
, db4
, db5
, db6
, db7
, db8
, db9
;
234 bAtStart
= (ftell(fp
) == 0);
236 clear_trxframe(fr
, FALSE
);
248 while (!fr
->bTitle
&& fgets2(line
, STRLEN
, fp
))
250 fr
->bTitle
= (strcmp(line
, "TITLE") == 0);
252 if (fr
->title
== NULL
)
254 fgets2(line
, STRLEN
, fp
);
255 fr
->title
= strdup(line
);
258 while (!bEnd
&& fgets2(line
, STRLEN
, fp
))
260 bEnd
= (strcmp(line
, "END") == 0);
262 fgets2(line
, STRLEN
, fp
);
265 /* Do not get a line if we are not at the start of the file, *
266 * because without a parameter file we don't know what is in *
267 * the trajectory and we have already read the line in the *
268 * previous call (VERY DIRTY). */
272 bTime
= (strcmp(line
, "TIMESTEP") == 0);
273 bAtoms
= (strcmp(line
, "POSITION") == 0);
274 bPos
= (bAtoms
|| (strcmp(line
, "POSITIONRED") == 0));
275 bVel
= (strncmp(line
, "VELOCITY", 8) == 0);
276 bBox
= (strcmp(line
, "BOX") == 0);
279 if (!fr
->bTime
&& !fr
->bX
)
285 bFinished
= (fgets2(line
, STRLEN
, fp
) == NULL
);
287 while (!bFinished
&& (line
[0] == '#'));
288 sscanf(line
, "%15d%15lf", &(fr
->step
), &db1
);
302 natoms
= read_g96_pos(line
, symtab
, fp
, infile
, fr
);
312 natoms
= read_g96_vel(line
, fp
, infile
, fr
);
319 while (!bEnd
&& fgets2(line
, STRLEN
, fp
))
321 bEnd
= (strncmp(line
, "END", 3) == 0);
322 if (!bEnd
&& (line
[0] != '#'))
324 nbp
= sscanf(line
, "%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf",
325 &db1
, &db2
, &db3
, &db4
, &db5
, &db6
, &db7
, &db8
, &db9
);
328 gmx_fatal(FARGS
, "Found a BOX line, but no box in %s", infile
);
330 fr
->box
[XX
][XX
] = db1
;
331 fr
->box
[YY
][YY
] = db2
;
332 fr
->box
[ZZ
][ZZ
] = db3
;
335 fr
->box
[XX
][YY
] = db4
;
336 fr
->box
[XX
][ZZ
] = db5
;
337 fr
->box
[YY
][XX
] = db6
;
338 fr
->box
[YY
][ZZ
] = db7
;
339 fr
->box
[ZZ
][XX
] = db8
;
340 fr
->box
[ZZ
][YY
] = db9
;
347 while (!bFinished
&& fgets2(line
, STRLEN
, fp
));
356 void write_g96_conf(FILE *out
, t_trxframe
*fr
,
357 int nindex
, const atom_id
*index
)
375 fprintf(out
, "TITLE\n%s\nEND\n", fr
->title
);
377 if (fr
->bStep
|| fr
->bTime
)
379 /* Officially the time format is %15.9, which is not enough for 10 ns */
380 fprintf(out
, "TIMESTEP\n%15d%15.6f\nEND\n", fr
->step
, fr
->time
);
386 fprintf(out
, "POSITION\n");
387 for (i
= 0; i
< nout
; i
++)
397 fprintf(out
, "%5d %-5s %-5s%7d%15.9f%15.9f%15.9f\n",
398 (atoms
->resinfo
[atoms
->atom
[a
].resind
].nr
) % 100000,
399 *atoms
->resinfo
[atoms
->atom
[a
].resind
].name
,
400 *atoms
->atomname
[a
], (i
+1) % 10000000,
401 fr
->x
[a
][XX
], fr
->x
[a
][YY
], fr
->x
[a
][ZZ
]);
406 fprintf(out
, "POSITIONRED\n");
407 for (i
= 0; i
< nout
; i
++)
417 fprintf(out
, "%15.9f%15.9f%15.9f\n",
418 fr
->x
[a
][XX
], fr
->x
[a
][YY
], fr
->x
[a
][ZZ
]);
421 fprintf(out
, "END\n");
427 fprintf(out
, "VELOCITY\n");
428 for (i
= 0; i
< nout
; i
++)
438 fprintf(out
, "%5d %-5s %-5s%7d%15.9f%15.9f%15.9f\n",
439 (atoms
->resinfo
[atoms
->atom
[a
].resind
].nr
) % 100000,
440 *atoms
->resinfo
[atoms
->atom
[a
].resind
].name
,
441 *atoms
->atomname
[a
], (i
+1) % 10000000,
442 fr
->v
[a
][XX
], fr
->v
[a
][YY
], fr
->v
[a
][ZZ
]);
447 fprintf(out
, "VELOCITYRED\n");
448 for (i
= 0; i
< nout
; i
++)
458 fprintf(out
, "%15.9f%15.9f%15.9f\n",
459 fr
->v
[a
][XX
], fr
->v
[a
][YY
], fr
->v
[a
][ZZ
]);
462 fprintf(out
, "END\n");
466 fprintf(out
, "BOX\n");
467 fprintf(out
, "%15.9f%15.9f%15.9f",
468 fr
->box
[XX
][XX
], fr
->box
[YY
][YY
], fr
->box
[ZZ
][ZZ
]);
469 if (fr
->box
[XX
][YY
] || fr
->box
[XX
][ZZ
] || fr
->box
[YY
][XX
] ||
470 fr
->box
[YY
][ZZ
] || fr
->box
[ZZ
][XX
] || fr
->box
[ZZ
][YY
])
472 fprintf(out
, "%15.9f%15.9f%15.9f%15.9f%15.9f%15.9f",
473 fr
->box
[XX
][YY
], fr
->box
[XX
][ZZ
], fr
->box
[YY
][XX
],
474 fr
->box
[YY
][ZZ
], fr
->box
[ZZ
][XX
], fr
->box
[ZZ
][YY
]);
477 fprintf(out
, "END\n");
481 static int get_espresso_word(FILE *fp
, char word
[])
493 if (i
== ' ' || i
== '\n' || i
== '\t')
518 word
[nc
++] = (char)i
;
522 while (i
!= EOF
&& ret
== 0);
526 /* printf("'%s'\n",word); */
531 static int check_open_parenthesis(FILE *fp
, int r
,
532 const char *infile
, const char *keyword
)
544 r
= get_espresso_word(fp
, word
);
551 gmx_fatal(FARGS
, "Expected '{' after '%s' in file '%s'",
559 static int check_close_parenthesis(FILE *fp
, int r
,
560 const char *infile
, const char *keyword
)
572 r
= get_espresso_word(fp
, word
);
579 gmx_fatal(FARGS
, "Expected '}' after section '%s' in file '%s'",
588 espID
, espPOS
, espTYPE
, espQ
, espV
, espF
, espMOLECULE
, espNR
590 const char *esp_prop
[espNR
] = {
591 "id", "pos", "type", "q", "v", "f",
595 static void read_espresso_conf(const char *infile
,
596 t_atoms
*atoms
, rvec x
[], rvec
*v
, matrix box
)
598 t_symtab
*symtab
= NULL
;
600 char word
[STRLEN
], buf
[STRLEN
];
601 int natoms
, level
, npar
, r
, nprop
, p
, i
, m
, molnr
;
604 gmx_bool bFoundParticles
, bFoundProp
, bFoundVariable
, bMol
;
614 fp
= gmx_fio_fopen(infile
, "r");
616 bFoundParticles
= FALSE
;
617 bFoundVariable
= FALSE
;
620 while ((r
= get_espresso_word(fp
, word
)))
622 if (level
== 1 && strcmp(word
, "particles") == 0 && !bFoundParticles
)
624 bFoundParticles
= TRUE
;
625 level
+= check_open_parenthesis(fp
, r
, infile
, "particles");
627 while (level
== 2 && (r
= get_espresso_word(fp
, word
)))
630 for (p
= 0; p
< espNR
; p
++)
632 if (strcmp(word
, esp_prop
[p
]) == 0)
636 /* printf(" prop[%d] = %s\n",nprop-1,esp_prop[prop[nprop-1]]); */
639 if (!bFoundProp
&& word
[0] != '}')
641 gmx_fatal(FARGS
, "Can not read Espresso files with particle property '%s'", word
);
643 if (bFoundProp
&& p
== espMOLECULE
)
654 while (level
> 0 && (r
= get_espresso_word(fp
, word
)))
666 for (p
= 0; p
< nprop
; p
++)
671 r
= get_espresso_word(fp
, word
);
675 for (m
= 0; m
< 3; m
++)
677 r
= get_espresso_word(fp
, word
);
678 sscanf(word
, "%lf", &d
);
683 r
= get_espresso_word(fp
, word
);
684 atoms
->atom
[i
].type
= strtol(word
, NULL
, 10);
687 r
= get_espresso_word(fp
, word
);
688 sscanf(word
, "%lf", &d
);
689 atoms
->atom
[i
].q
= d
;
692 for (m
= 0; m
< 3; m
++)
694 r
= get_espresso_word(fp
, word
);
695 sscanf(word
, "%lf", &d
);
700 for (m
= 0; m
< 3; m
++)
702 r
= get_espresso_word(fp
, word
);
707 r
= get_espresso_word(fp
, word
);
708 molnr
= strtol(word
, NULL
, 10);
710 atoms
->resinfo
[atoms
->atom
[i
-1].resind
].nr
!= molnr
)
712 atoms
->atom
[i
].resind
=
713 (i
== 0 ? 0 : atoms
->atom
[i
-1].resind
+1);
714 atoms
->resinfo
[atoms
->atom
[i
].resind
].nr
= molnr
;
715 atoms
->resinfo
[atoms
->atom
[i
].resind
].ic
= ' ';
716 atoms
->resinfo
[atoms
->atom
[i
].resind
].chainid
= ' ';
717 atoms
->resinfo
[atoms
->atom
[i
].resind
].chainnum
= molnr
; /* Not sure if this is right? */
721 atoms
->atom
[i
].resind
= atoms
->atom
[i
-1].resind
;
726 /* Generate an atom name from the particle type */
727 sprintf(buf
, "T%d", atoms
->atom
[i
].type
);
728 atoms
->atomname
[i
] = put_symtab(symtab
, buf
);
731 if (i
== 0 || atoms
->atom
[i
].resind
!= atoms
->atom
[i
-1].resind
)
733 atoms
->resinfo
[atoms
->atom
[i
].resind
].name
=
734 put_symtab(symtab
, "MOL");
739 /* Residue number is the atom number */
740 atoms
->atom
[i
].resind
= i
;
741 /* Generate an residue name from the particle type */
742 if (atoms
->atom
[i
].type
< 26)
744 sprintf(buf
, "T%c", 'A'+atoms
->atom
[i
].type
);
748 sprintf(buf
, "T%c%c",
749 'A'+atoms
->atom
[i
].type
/26, 'A'+atoms
->atom
[i
].type
%26);
751 t_atoms_set_resinfo(atoms
, i
, symtab
, buf
, i
, ' ', 0, ' ');
761 atoms
->nres
= atoms
->nr
;
765 gmx_fatal(FARGS
, "Internal inconsistency in Espresso routines, read %d atoms, expected %d atoms", i
, atoms
->nr
);
768 else if (level
== 1 && strcmp(word
, "variable") == 0 && !bFoundVariable
)
770 bFoundVariable
= TRUE
;
771 level
+= check_open_parenthesis(fp
, r
, infile
, "variable");
772 while (level
== 2 && (r
= get_espresso_word(fp
, word
)))
774 if (level
== 2 && strcmp(word
, "box_l") == 0)
776 for (m
= 0; m
< 3; m
++)
778 r
= get_espresso_word(fp
, word
);
779 sscanf(word
, "%lf", &d
);
782 level
+= check_close_parenthesis(fp
, r
, infile
, "box_l");
796 if (!bFoundParticles
)
798 fprintf(stderr
, "Did not find a particles section in Espresso file '%s'\n",
805 static int get_espresso_coordnum(const char *infile
)
809 int natoms
, level
, r
;
810 gmx_bool bFoundParticles
;
814 fp
= gmx_fio_fopen(infile
, "r");
816 bFoundParticles
= FALSE
;
818 while ((r
= get_espresso_word(fp
, word
)) && !bFoundParticles
)
820 if (level
== 1 && strcmp(word
, "particles") == 0 && !bFoundParticles
)
822 bFoundParticles
= TRUE
;
823 level
+= check_open_parenthesis(fp
, r
, infile
, "particles");
824 while (level
> 0 && (r
= get_espresso_word(fp
, word
)))
849 if (!bFoundParticles
)
851 fprintf(stderr
, "Did not find a particles section in Espresso file '%s'\n",
860 static void write_espresso_conf_indexed(FILE *out
, const char *title
,
861 t_atoms
*atoms
, int nx
, atom_id
*index
,
862 rvec
*x
, rvec
*v
, matrix box
)
866 fprintf(out
, "# %s\n", title
);
869 gmx_warning("The Espresso format does not support triclinic unit-cells");
871 fprintf(out
, "{variable {box_l %f %f %f}}\n", box
[0][0], box
[1][1], box
[2][2]);
873 fprintf(out
, "{particles {id pos type q%s}\n", v
? " v" : "");
874 for (i
= 0; i
< nx
; i
++)
884 fprintf(out
, "\t{%d %f %f %f %d %g",
885 j
, x
[j
][XX
], x
[j
][YY
], x
[j
][ZZ
],
886 atoms
->atom
[j
].type
, atoms
->atom
[j
].q
);
889 fprintf(out
, " %f %f %f", v
[j
][XX
], v
[j
][YY
], v
[j
][ZZ
]);
896 static void get_coordnum_fp (FILE *in
, char *title
, int *natoms
)
900 fgets2 (title
, STRLEN
, in
);
901 fgets2 (line
, STRLEN
, in
);
902 if (sscanf (line
, "%d", natoms
) != 1)
904 gmx_fatal(FARGS
, "gro file does not have the number of atoms on the second line");
908 static void get_coordnum (const char *infile
, int *natoms
)
913 in
= gmx_fio_fopen(infile
, "r");
914 get_coordnum_fp(in
, title
, natoms
);
918 static gmx_bool
get_w_conf(FILE *in
, const char *infile
, char *title
,
919 t_symtab
*symtab
, t_atoms
*atoms
, int *ndec
,
920 rvec x
[], rvec
*v
, matrix box
)
923 char line
[STRLEN
+1], *ptr
;
925 double x1
, y1
, z1
, x2
, y2
, z2
;
927 int natoms
, i
, m
, resnr
, newres
, oldres
, ddist
, c
;
928 gmx_bool bFirst
, bVel
;
932 oldres
= NOTSET
; /* Unlikely number for the first residue! */
935 /* Read the title and number of atoms */
936 get_coordnum_fp(in
, title
, &natoms
);
938 if (natoms
> atoms
->nr
)
940 gmx_fatal(FARGS
, "gro file contains more atoms (%d) than expected (%d)",
943 else if (natoms
< atoms
->nr
)
945 fprintf(stderr
, "Warning: gro file contains less atoms (%d) than expected"
946 " (%d)\n", natoms
, atoms
->nr
);
953 /* just pray the arrays are big enough */
954 for (i
= 0; (i
< natoms
); i
++)
956 if ((fgets2 (line
, STRLEN
, in
)) == NULL
)
958 gmx_fatal(FARGS
, "Unexpected end of file in file %s at line %d",
961 if (strlen(line
) < 39)
963 gmx_fatal(FARGS
, "Invalid line in %s for atom %d:\n%s", infile
, i
+1, line
);
966 /* determine read precision from distance between periods
971 p1
= strchr(line
, '.');
974 gmx_fatal(FARGS
, "A coordinate in file %s does not contain a '.'", infile
);
976 p2
= strchr(&p1
[1], '.');
979 gmx_fatal(FARGS
, "A coordinate in file %s does not contain a '.'", infile
);
984 p3
= strchr(&p2
[1], '.');
987 gmx_fatal(FARGS
, "A coordinate in file %s does not contain a '.'", infile
);
990 if (p3
- p2
!= ddist
)
992 gmx_fatal(FARGS
, "The spacing of the decimal points in file %s is not consistent for x, y and z", infile
);
997 memcpy(name
, line
, 5);
999 sscanf(name
, "%d", &resnr
);
1000 memcpy(name
, line
+5, 5);
1002 if (resnr
!= oldres
)
1006 if (newres
>= natoms
)
1008 gmx_fatal(FARGS
, "More residues than atoms in %s (natoms = %d)",
1011 atoms
->atom
[i
].resind
= newres
;
1012 t_atoms_set_resinfo(atoms
, i
, symtab
, name
, resnr
, ' ', 0, ' ');
1016 atoms
->atom
[i
].resind
= newres
;
1020 memcpy(name
, line
+10, 5);
1021 atoms
->atomname
[i
] = put_symtab(symtab
, name
);
1023 /* eventueel controle atomnumber met i+1 */
1025 /* coordinates (start after residue data) */
1027 /* Read fixed format */
1028 for (m
= 0; m
< DIM
; m
++)
1030 for (c
= 0; (c
< ddist
&& ptr
[0]); c
++)
1036 if (sscanf (buf
, "%lf %lf", &x1
, &x2
) != 1)
1038 gmx_fatal(FARGS
, "Something is wrong in the coordinate formatting of file %s. Note that gro is fixed format (see the manual)", infile
);
1046 /* velocities (start after residues and coordinates) */
1049 /* Read fixed format */
1050 for (m
= 0; m
< DIM
; m
++)
1052 for (c
= 0; (c
< ddist
&& ptr
[0]); c
++)
1058 if (sscanf (buf
, "%lf", &x1
) != 1)
1070 atoms
->nres
= newres
+ 1;
1073 fgets2 (line
, STRLEN
, in
);
1074 if (sscanf (line
, "%lf%lf%lf", &x1
, &y1
, &z1
) != 3)
1076 gmx_warning("Bad box in file %s", infile
);
1078 /* Generate a cubic box */
1079 for (m
= 0; (m
< DIM
); m
++)
1081 xmin
[m
] = xmax
[m
] = x
[0][m
];
1083 for (i
= 1; (i
< atoms
->nr
); i
++)
1085 for (m
= 0; (m
< DIM
); m
++)
1087 xmin
[m
] = min(xmin
[m
], x
[i
][m
]);
1088 xmax
[m
] = max(xmax
[m
], x
[i
][m
]);
1091 for (i
= 0; i
< DIM
; i
++)
1093 for (m
= 0; m
< DIM
; m
++)
1098 for (m
= 0; (m
< DIM
); m
++)
1100 box
[m
][m
] = (xmax
[m
]-xmin
[m
]);
1102 fprintf(stderr
, "Generated a cubic box %8.3f x %8.3f x %8.3f\n",
1103 box
[XX
][XX
], box
[YY
][YY
], box
[ZZ
][ZZ
]);
1107 /* We found the first three values, the diagonal elements */
1111 if (sscanf (line
, "%*f%*f%*f%lf%lf%lf%lf%lf%lf",
1112 &x1
, &y1
, &z1
, &x2
, &y2
, &z2
) != 6)
1114 x1
= y1
= z1
= x2
= y2
= z2
= 0.0;
1127 static void read_whole_conf(const char *infile
, char *title
,
1128 t_atoms
*atoms
, rvec x
[], rvec
*v
, matrix box
)
1135 in
= gmx_fio_fopen(infile
, "r");
1137 open_symtab(&symtab
);
1138 get_w_conf(in
, infile
, title
, &symtab
, atoms
, &ndec
, x
, v
, box
);
1139 /* We can't free the symbols, as they are still used in atoms, so
1140 * the only choice is to leak them. */
1141 free_symtab(&symtab
);
1146 static gmx_bool
gmx_one_before_eof(FILE *fp
)
1151 if ((beof
= fread(data
, 1, 1, fp
)) == 1)
1153 gmx_fseek(fp
, -1, SEEK_CUR
);
1158 gmx_bool
gro_next_x_or_v(FILE *status
, t_trxframe
*fr
)
1162 char title
[STRLEN
], *p
;
1166 if (gmx_one_before_eof(status
))
1171 open_symtab(&symtab
);
1172 atoms
.nr
= fr
->natoms
;
1173 snew(atoms
.atom
, fr
->natoms
);
1174 atoms
.nres
= fr
->natoms
;
1175 snew(atoms
.resinfo
, fr
->natoms
);
1176 snew(atoms
.atomname
, fr
->natoms
);
1178 fr
->bV
= get_w_conf(status
, title
, title
, &symtab
, &atoms
, &ndec
, fr
->x
, fr
->v
, fr
->box
);
1181 /* prec = 10^ndec: */
1182 for (i
= 0; i
< ndec
; i
++)
1192 sfree(atoms
.resinfo
);
1193 sfree(atoms
.atomname
);
1194 done_symtab(&symtab
);
1196 if ((p
= strstr(title
, "t=")) != NULL
)
1199 if (sscanf(p
, "%lf", &tt
) == 1)
1211 if (atoms
.nr
!= fr
->natoms
)
1213 gmx_fatal(FARGS
, "Number of atoms in gro frame (%d) doesn't match the number in the previous frame (%d)", atoms
.nr
, fr
->natoms
);
1219 int gro_first_x_or_v(FILE *status
, t_trxframe
*fr
)
1224 fprintf(stderr
, "Reading frames from gro file");
1225 get_coordnum_fp(status
, title
, &fr
->natoms
);
1227 fprintf(stderr
, " '%s', %d atoms.\n", title
, fr
->natoms
);
1230 if (fr
->natoms
== 0)
1232 gmx_file("No coordinates in gro file");
1235 snew(fr
->x
, fr
->natoms
);
1236 snew(fr
->v
, fr
->natoms
);
1237 gro_next_x_or_v(status
, fr
);
1242 static void make_hconf_format(int pr
, gmx_bool bVel
, char format
[])
1246 /* build format string for printing,
1247 something like "%8.3f" for x and "%8.4f" for v */
1260 sprintf(format
, "%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df\n",
1261 l
, pr
, l
, pr
, l
, pr
, l
, vpr
, l
, vpr
, l
, vpr
);
1265 sprintf(format
, "%%%d.%df%%%d.%df%%%d.%df\n", l
, pr
, l
, pr
, l
, pr
);
1270 static void write_hconf_box(FILE *out
, int pr
, matrix box
)
1281 if (box
[XX
][YY
] || box
[XX
][ZZ
] || box
[YY
][XX
] || box
[YY
][ZZ
] ||
1282 box
[ZZ
][XX
] || box
[ZZ
][YY
])
1284 sprintf(format
, "%%%d.%df%%%d.%df%%%d.%df"
1285 "%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df\n",
1286 l
, pr
, l
, pr
, l
, pr
, l
, pr
, l
, pr
, l
, pr
, l
, pr
, l
, pr
, l
, pr
);
1287 fprintf(out
, format
,
1288 box
[XX
][XX
], box
[YY
][YY
], box
[ZZ
][ZZ
],
1289 box
[XX
][YY
], box
[XX
][ZZ
], box
[YY
][XX
],
1290 box
[YY
][ZZ
], box
[ZZ
][XX
], box
[ZZ
][YY
]);
1294 sprintf(format
, "%%%d.%df%%%d.%df%%%d.%df\n", l
, pr
, l
, pr
, l
, pr
);
1295 fprintf(out
, format
,
1296 box
[XX
][XX
], box
[YY
][YY
], box
[ZZ
][ZZ
]);
1300 void write_hconf_indexed_p(FILE *out
, const char *title
, t_atoms
*atoms
,
1301 int nx
, const atom_id index
[], int pr
,
1302 rvec
*x
, rvec
*v
, matrix box
)
1304 char resnm
[6], nm
[6], format
[100];
1305 int ai
, i
, resind
, resnr
;
1307 bromacs(format
, 99);
1308 fprintf (out
, "%s\n", (title
&& title
[0]) ? title
: format
);
1309 fprintf (out
, "%5d\n", nx
);
1311 make_hconf_format(pr
, v
!= NULL
, format
);
1313 for (i
= 0; (i
< nx
); i
++)
1317 resind
= atoms
->atom
[ai
].resind
;
1318 strncpy(resnm
, " ??? ", sizeof(resnm
)-1);
1319 if (resind
< atoms
->nres
)
1321 strncpy(resnm
, *atoms
->resinfo
[resind
].name
, sizeof(resnm
)-1);
1322 resnr
= atoms
->resinfo
[resind
].nr
;
1326 strncpy(resnm
, " ??? ", sizeof(resnm
)-1);
1332 strncpy(nm
, *atoms
->atomname
[ai
], sizeof(nm
)-1);
1336 strncpy(nm
, " ??? ", sizeof(nm
)-1);
1339 fprintf(out
, "%5d%-5.5s%5.5s%5d", resnr
%100000, resnm
, nm
, (ai
+1)%100000);
1340 /* next fprintf uses built format string */
1343 fprintf(out
, format
,
1344 x
[ai
][XX
], x
[ai
][YY
], x
[ai
][ZZ
], v
[ai
][XX
], v
[ai
][YY
], v
[ai
][ZZ
]);
1348 fprintf(out
, format
,
1349 x
[ai
][XX
], x
[ai
][YY
], x
[ai
][ZZ
]);
1353 write_hconf_box(out
, pr
, box
);
1358 static void write_hconf_mtop(FILE *out
, const char *title
, gmx_mtop_t
*mtop
,
1360 rvec
*x
, rvec
*v
, matrix box
)
1364 gmx_mtop_atomloop_all_t aloop
;
1366 char *atomname
, *resname
;
1368 bromacs(format
, 99);
1369 fprintf (out
, "%s\n", (title
&& title
[0]) ? title
: format
);
1370 fprintf (out
, "%5d\n", mtop
->natoms
);
1372 make_hconf_format(pr
, v
!= NULL
, format
);
1374 aloop
= gmx_mtop_atomloop_all_init(mtop
);
1375 while (gmx_mtop_atomloop_all_next(aloop
, &i
, &atom
))
1377 gmx_mtop_atomloop_all_names(aloop
, &atomname
, &resnr
, &resname
);
1379 fprintf(out
, "%5d%-5.5s%5.5s%5d",
1380 resnr
%100000, resname
, atomname
, (i
+1)%100000);
1381 /* next fprintf uses built format string */
1384 fprintf(out
, format
,
1385 x
[i
][XX
], x
[i
][YY
], x
[i
][ZZ
], v
[i
][XX
], v
[i
][YY
], v
[i
][ZZ
]);
1389 fprintf(out
, format
,
1390 x
[i
][XX
], x
[i
][YY
], x
[i
][ZZ
]);
1394 write_hconf_box(out
, pr
, box
);
1399 void write_hconf_p(FILE *out
, const char *title
, t_atoms
*atoms
, int pr
,
1400 rvec
*x
, rvec
*v
, matrix box
)
1405 snew(aa
, atoms
->nr
);
1406 for (i
= 0; (i
< atoms
->nr
); i
++)
1410 write_hconf_indexed_p(out
, title
, atoms
, atoms
->nr
, aa
, pr
, x
, v
, box
);
1414 void write_conf_p(const char *outfile
, const char *title
,
1415 t_atoms
*atoms
, int pr
,
1416 rvec
*x
, rvec
*v
, matrix box
)
1420 out
= gmx_fio_fopen(outfile
, "w");
1421 write_hconf_p(out
, title
, atoms
, pr
, x
, v
, box
);
1423 gmx_fio_fclose (out
);
1426 static void write_conf(const char *outfile
, const char *title
, t_atoms
*atoms
,
1427 rvec
*x
, rvec
*v
, matrix box
)
1429 write_conf_p(outfile
, title
, atoms
, 3, x
, v
, box
);
1432 void write_sto_conf_indexed(const char *outfile
, const char *title
,
1434 rvec x
[], rvec
*v
, int ePBC
, matrix box
,
1435 atom_id nindex
, atom_id index
[])
1441 ftp
= fn2ftp(outfile
);
1445 out
= gmx_fio_fopen(outfile
, "w");
1446 write_hconf_indexed_p(out
, title
, atoms
, nindex
, index
, 3, x
, v
, box
);
1447 gmx_fio_fclose(out
);
1450 clear_trxframe(&fr
, TRUE
);
1453 fr
.natoms
= atoms
->nr
;
1464 copy_mat(box
, fr
.box
);
1465 out
= gmx_fio_fopen(outfile
, "w");
1466 write_g96_conf(out
, &fr
, nindex
, index
);
1467 gmx_fio_fclose(out
);
1473 out
= gmx_fio_fopen(outfile
, "w");
1474 write_pdbfile_indexed(out
, title
, atoms
, x
, ePBC
, box
, ' ', -1, nindex
, index
, NULL
, TRUE
);
1475 gmx_fio_fclose(out
);
1478 out
= gmx_fio_fopen(outfile
, "w");
1479 write_espresso_conf_indexed(out
, title
, atoms
, nindex
, index
, x
, v
, box
);
1480 gmx_fio_fclose(out
);
1485 gmx_fatal(FARGS
, "Sorry, can not write a topology to %s", outfile
);
1488 gmx_incons("Not supported in write_sto_conf_indexed");
1492 void write_sto_conf(const char *outfile
, const char *title
, t_atoms
*atoms
,
1493 rvec x
[], rvec
*v
, int ePBC
, matrix box
)
1499 ftp
= fn2ftp(outfile
);
1503 write_conf(outfile
, title
, atoms
, x
, v
, box
);
1506 clear_trxframe(&fr
, TRUE
);
1509 fr
.natoms
= atoms
->nr
;
1520 copy_mat(box
, fr
.box
);
1521 out
= gmx_fio_fopen(outfile
, "w");
1522 write_g96_conf(out
, &fr
, -1, NULL
);
1523 gmx_fio_fclose(out
);
1528 out
= gmx_fio_fopen(outfile
, "w");
1529 write_pdbfile(out
, title
, atoms
, x
, ePBC
, box
, ' ', -1, NULL
, TRUE
);
1530 gmx_fio_fclose(out
);
1533 out
= gmx_fio_fopen(outfile
, "w");
1534 write_espresso_conf_indexed(out
, title
, atoms
, atoms
->nr
, NULL
, x
, v
, box
);
1535 gmx_fio_fclose(out
);
1540 gmx_fatal(FARGS
, "Sorry, can not write a topology to %s", outfile
);
1543 gmx_incons("Not supported in write_sto_conf");
1547 void write_sto_conf_mtop(const char *outfile
, const char *title
,
1549 rvec x
[], rvec
*v
, int ePBC
, matrix box
)
1555 ftp
= fn2ftp(outfile
);
1559 out
= gmx_fio_fopen(outfile
, "w");
1560 write_hconf_mtop(out
, title
, mtop
, 3, x
, v
, box
);
1561 gmx_fio_fclose(out
);
1564 /* This is a brute force approach which requires a lot of memory.
1565 * We should implement mtop versions of all writing routines.
1567 atoms
= gmx_mtop_global_atoms(mtop
);
1569 write_sto_conf(outfile
, title
, &atoms
, x
, v
, ePBC
, box
);
1576 void get_stx_coordnum(const char *infile
, int *natoms
)
1579 int ftp
, tpxver
, tpxgen
;
1581 char g96_line
[STRLEN
+1];
1583 ftp
= fn2ftp(infile
);
1584 range_check(ftp
, 0, efNR
);
1588 get_coordnum(infile
, natoms
);
1591 in
= gmx_fio_fopen(infile
, "r");
1598 *natoms
= read_g96_conf(in
, infile
, &fr
, g96_line
);
1604 in
= gmx_fio_fopen(infile
, "r");
1605 get_pdb_coordnum(in
, natoms
);
1609 *natoms
= get_espresso_coordnum(infile
);
1617 read_tpxheader(infile
, &tpx
, TRUE
, &tpxver
, &tpxgen
);
1618 *natoms
= tpx
.natoms
;
1622 gmx_fatal(FARGS
, "File type %s not supported in get_stx_coordnum",
1627 void read_stx_conf(const char *infile
, char *title
, t_atoms
*atoms
,
1628 rvec x
[], rvec
*v
, int *ePBC
, matrix box
)
1637 char g96_line
[STRLEN
+1];
1641 fprintf(stderr
, "Warning: Number of atoms in %s is 0\n", infile
);
1643 else if (atoms
->atom
== NULL
)
1645 gmx_mem("Uninitialized array atom");
1653 ftp
= fn2ftp(infile
);
1657 read_whole_conf(infile
, title
, atoms
, x
, v
, box
);
1661 fr
.natoms
= atoms
->nr
;
1666 in
= gmx_fio_fopen(infile
, "r");
1667 read_g96_conf(in
, infile
, &fr
, g96_line
);
1669 copy_mat(fr
.box
, box
);
1670 strncpy(title
, fr
.title
, STRLEN
);
1675 read_pdb_conf(infile
, title
, atoms
, x
, ePBC
, box
, TRUE
, NULL
);
1678 read_espresso_conf(infile
, atoms
, x
, v
, box
);
1684 i
= read_tpx(infile
, NULL
, box
, &natoms
, x
, v
, NULL
, mtop
);
1690 strcpy(title
, *(mtop
->name
));
1692 /* Free possibly allocated memory */
1695 *atoms
= gmx_mtop_global_atoms(mtop
);
1696 top
= gmx_mtop_t_to_t_topology(mtop
);
1697 tpx_make_chain_identifiers(atoms
, &top
.mols
);
1700 /* The strings in the symtab are still in use in the returned t_atoms
1701 * structure, so we should not free them. But there is no place to put the
1702 * symbols; the only choice is to leak the memory...
1703 * So we clear the symbol table before freeing the topology structure. */
1704 free_symtab(&top
.symtab
);
1709 gmx_incons("Not supported in read_stx_conf");