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
60 typedef struct gmx_conect_t {
63 gmx_conection_t *conect;
66 static const char *pdbtp[epdbNR] = {
67 "ATOM ", "HETATM", "ANISOU", "CRYST1",
68 "COMPND", "MODEL", "ENDMDL", "TER", "HEADER", "TITLE", "REMARK",
73 /* this is not very good,
74 but these are only used in gmx_trjconv and gmx_editconv */
75 static gmx_bool bWideFormat = FALSE;
76 #define REMARK_SIM_BOX "REMARK THIS IS A SIMULATION BOX"
78 void set_pdb_wide_format(gmx_bool bSet)
83 static void xlate_atomname_pdb2gmx(char *name)
88 length = strlen(name);
89 if (length > 3 && isdigit(name[0]))
92 for (i = 1; i < length; i++)
96 name[length-1] = temp;
100 static void xlate_atomname_gmx2pdb(char *name)
105 length = strlen(name);
106 if (length > 3 && isdigit(name[length-1]))
108 temp = name[length-1];
109 for (i = length-1; i > 0; --i)
118 void gmx_write_pdb_box(FILE *out, int ePBC, matrix box)
120 real alpha, beta, gamma;
124 ePBC = guess_ePBC(box);
127 if (ePBC == epbcNONE)
132 if (norm2(box[YY])*norm2(box[ZZ]) != 0)
134 alpha = RAD2DEG*acos(cos_angle_no_table(box[YY], box[ZZ]));
140 if (norm2(box[XX])*norm2(box[ZZ]) != 0)
142 beta = RAD2DEG*acos(cos_angle_no_table(box[XX], box[ZZ]));
148 if (norm2(box[XX])*norm2(box[YY]) != 0)
150 gamma = RAD2DEG*acos(cos_angle_no_table(box[XX], box[YY]));
156 fprintf(out, "REMARK THIS IS A SIMULATION BOX\n");
157 if (ePBC != epbcSCREW)
159 fprintf(out, "CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f %-11s%4d\n",
160 10*norm(box[XX]), 10*norm(box[YY]), 10*norm(box[ZZ]),
161 alpha, beta, gamma, "P 1", 1);
165 /* Double the a-vector length and write the correct space group */
166 fprintf(out, "CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f %-11s%4d\n",
167 20*norm(box[XX]), 10*norm(box[YY]), 10*norm(box[ZZ]),
168 alpha, beta, gamma, "P 21 1 1", 1);
173 static void read_cryst1(char *line, int *ePBC, matrix box)
176 char sa[12], sb[12], sc[12], sg[SG_SIZE+1], ident;
177 double fa, fb, fc, alpha, beta, gamma, cosa, cosb, cosg, sing;
178 int syma, symb, symc;
181 sscanf(line, "%*s%s%s%s%lf%lf%lf", sa, sb, sc, &alpha, &beta, &gamma);
184 if (strlen(line) >= 55)
186 strncpy(sg, line+55, SG_SIZE);
192 sscanf(sg, "%c %d %d %d", &ident, &syma, &symb, &symc);
193 if (ident == 'P' && syma == 1 && symb <= 1 && symc <= 1)
195 fc = strtod(sc, NULL)*0.1;
196 ePBC_file = (fc > 0 ? epbcXYZ : epbcXY);
198 if (ident == 'P' && syma == 21 && symb == 1 && symc == 1)
200 ePBC_file = epbcSCREW;
210 fa = strtod(sa, NULL)*0.1;
211 fb = strtod(sb, NULL)*0.1;
212 fc = strtod(sc, NULL)*0.1;
213 if (ePBC_file == epbcSCREW)
219 if ((alpha != 90.0) || (beta != 90.0) || (gamma != 90.0))
223 cosa = cos(alpha*DEG2RAD);
231 cosb = cos(beta*DEG2RAD);
239 cosg = cos(gamma*DEG2RAD);
240 sing = sin(gamma*DEG2RAD);
247 box[YY][XX] = fb*cosg;
248 box[YY][YY] = fb*sing;
249 box[ZZ][XX] = fc*cosb;
250 box[ZZ][YY] = fc*(cosa - cosb*cosg)/sing;
251 box[ZZ][ZZ] = sqrt(fc*fc
252 - box[ZZ][XX]*box[ZZ][XX] - box[ZZ][YY]*box[ZZ][YY]);
262 void write_pdbfile_indexed(FILE *out, const char *title,
263 t_atoms *atoms, rvec x[],
264 int ePBC, matrix box, char chainid,
265 int model_nr, atom_id nindex, const atom_id index[],
266 gmx_conect conect, gmx_bool bTerSepChains)
268 gmx_conect_t *gc = (gmx_conect_t *)conect;
269 char resnm[6], nm[6], pdbform[128], pukestring[100];
271 int resind, resnr, type;
272 unsigned char resic, ch;
276 int chainnum, lastchainnum;
277 int lastresind, lastchainresind;
278 gmx_residuetype_t rt;
279 const char *p_restype;
280 const char *p_lastrestype;
282 gmx_residuetype_init(&rt);
284 bromacs(pukestring, 99);
285 fprintf(out, "TITLE %s\n", (title && title[0]) ? title : pukestring);
288 fprintf(out, "REMARK This file does not adhere to the PDB standard\n");
289 fprintf(out, "REMARK As a result of, some programs may not like it\n");
291 if (box && ( norm2(box[XX]) || norm2(box[YY]) || norm2(box[ZZ]) ) )
293 gmx_write_pdb_box(out, ePBC, box);
297 /* Check whether any occupancies are set, in that case leave it as is,
298 * otherwise set them all to one
301 for (ii = 0; (ii < nindex) && bOccup; ii++)
304 bOccup = bOccup && (atoms->pdbinfo[i].occup == 0.0);
312 fprintf(out, "MODEL %8d\n", model_nr > 0 ? model_nr : 1);
314 lastchainresind = -1;
319 for (ii = 0; ii < nindex; ii++)
323 resind = atoms->atom[i].resind;
324 chainnum = atoms->resinfo[resind].chainnum;
325 p_lastrestype = p_restype;
326 gmx_residuetype_get_type(rt, *atoms->resinfo[resind].name, &p_restype);
328 /* Add a TER record if we changed chain, and if either the previous or this chain is protein/DNA/RNA. */
329 if (bTerSepChains && ii > 0 && chainnum != lastchainnum)
331 /* Only add TER if the previous chain contained protein/DNA/RNA. */
332 if (gmx_residuetype_is_protein(rt, p_lastrestype) || gmx_residuetype_is_dna(rt, p_lastrestype) || gmx_residuetype_is_rna(rt, p_lastrestype))
334 fprintf(out, "TER\n");
336 lastchainnum = chainnum;
337 lastchainresind = lastresind;
340 strncpy(resnm, *atoms->resinfo[resind].name, sizeof(resnm)-1);
341 strncpy(nm, *atoms->atomname[i], sizeof(nm)-1);
342 /* rename HG12 to 2HG1, etc. */
343 xlate_atomname_gmx2pdb(nm);
344 resnr = atoms->resinfo[resind].nr;
345 resic = atoms->resinfo[resind].ic;
352 ch = atoms->resinfo[resind].chainid;
361 resnr = resnr % 10000;
365 type = atoms->pdbinfo[i].type;
366 occup = bOccup ? 1.0 : atoms->pdbinfo[i].occup;
367 bfac = atoms->pdbinfo[i].bfac;
378 "%-6s%5u %-4.4s %3.3s %c%4d%c %10.5f%10.5f%10.5f%8.4f%8.4f %2s\n");
382 /* Check whether atomname is an element name */
383 if ((strlen(nm) < 4) && (gmx_strcasecmp(nm, atoms->atom[i].elem) != 0))
385 strcpy(pdbform, get_pdbformat());
389 strcpy(pdbform, get_pdbformat4());
393 if (nlongname < maxwln)
395 fprintf(stderr, "WARNING: Writing out atom name (%s) longer than 4 characters to .pdb file\n", nm);
397 else if (nlongname == maxwln)
399 fprintf(stderr, "WARNING: More than %d long atom names, will not write more warnings\n", maxwln);
404 strcat(pdbform, "%6.2f%6.2f %2s\n");
406 fprintf(out, pdbform, pdbtp[type], (i+1)%100000, nm, resnm, ch, resnr,
407 (resic == '\0') ? ' ' : resic,
408 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], occup, bfac, atoms->atom[i].elem);
409 if (atoms->pdbinfo && atoms->pdbinfo[i].bAnisotropic)
411 fprintf(out, "ANISOU%5u %-4.4s%3.3s %c%4d%c %7d%7d%7d%7d%7d%7d\n",
412 (i+1)%100000, nm, resnm, ch, resnr,
413 (resic == '\0') ? ' ' : resic,
414 atoms->pdbinfo[i].uij[0], atoms->pdbinfo[i].uij[1],
415 atoms->pdbinfo[i].uij[2], atoms->pdbinfo[i].uij[3],
416 atoms->pdbinfo[i].uij[4], atoms->pdbinfo[i].uij[5]);
420 fprintf(out, "TER\n");
421 fprintf(out, "ENDMDL\n");
425 /* Write conect records */
426 for (i = 0; (i < gc->nconect); i++)
428 fprintf(out, "CONECT%5d%5d\n", gc->conect[i].ai+1, gc->conect[i].aj+1);
432 gmx_residuetype_destroy(rt);
435 void write_pdbfile(FILE *out, const char *title, t_atoms *atoms, rvec x[],
436 int ePBC, matrix box, char chainid, int model_nr, gmx_conect conect, gmx_bool bTerSepChains)
440 snew(index, atoms->nr);
441 for (i = 0; i < atoms->nr; i++)
445 write_pdbfile_indexed(out, title, atoms, x, ePBC, box, chainid, model_nr,
446 atoms->nr, index, conect, bTerSepChains);
450 int line2type(char *line)
455 for (k = 0; (k < 6); k++)
461 for (k = 0; (k < epdbNR); k++)
463 if (strncmp(type, pdbtp[k], strlen(pdbtp[k])) == 0)
472 static void read_anisou(char line[], int natom, t_atoms *atoms)
476 char anr[12], anm[12];
480 for (k = 0; (k < 5); k++, j++)
486 for (k = 0; (k < 4); k++, j++)
493 /* Strip off spaces */
496 /* Search backwards for number and name only */
497 atomnr = strtol(anr, NULL, 10);
498 for (i = natom-1; (i >= 0); i--)
500 if ((strcmp(anm, *(atoms->atomname[i])) == 0) &&
501 (atomnr == atoms->pdbinfo[i].atomnr))
508 fprintf(stderr, "Skipping ANISOU record (atom %s %d not found)\n",
513 if (sscanf(line+29, "%d%d%d%d%d%d",
514 &atoms->pdbinfo[i].uij[U11], &atoms->pdbinfo[i].uij[U22],
515 &atoms->pdbinfo[i].uij[U33], &atoms->pdbinfo[i].uij[U12],
516 &atoms->pdbinfo[i].uij[U13], &atoms->pdbinfo[i].uij[U23])
519 atoms->pdbinfo[i].bAnisotropic = TRUE;
523 fprintf(stderr, "Invalid ANISOU record for atom %d\n", i);
524 atoms->pdbinfo[i].bAnisotropic = FALSE;
529 void get_pdb_atomnumber(t_atoms *atoms, gmx_atomprop_t aps)
531 int i, atomnumber, len;
533 char anm[6], anm_copy[6], *ptr;
539 gmx_incons("Trying to deduce atomnumbers when no pdb information is present");
541 for (i = 0; (i < atoms->nr); i++)
543 strcpy(anm, atoms->pdbinfo[i].atomnm);
544 strcpy(anm_copy, atoms->pdbinfo[i].atomnm);
547 if ((anm[0] != ' ') && ((len <= 2) || ((len > 2) && !isdigit(anm[2]))))
550 if (gmx_atomprop_query(aps, epropElement, "???", anm_copy, &eval))
552 atomnumber = gmx_nint(eval);
557 if (gmx_atomprop_query(aps, epropElement, "???", anm_copy, &eval))
559 atomnumber = gmx_nint(eval);
563 if (atomnumber == NOTSET)
566 while ((k < strlen(anm)) && (isspace(anm[k]) || isdigit(anm[k])))
570 anm_copy[0] = anm[k];
572 if (gmx_atomprop_query(aps, epropElement, "???", anm_copy, &eval))
574 atomnumber = gmx_nint(eval);
577 atoms->atom[i].atomnumber = atomnumber;
578 ptr = gmx_atomprop_element(aps, atomnumber);
579 strncpy(atoms->atom[i].elem, ptr == NULL ? "" : ptr, 4);
582 fprintf(debug, "Atomnumber for atom '%s' is %d\n", anm, atomnumber);
587 static int read_atom(t_symtab *symtab,
588 char line[], int type, int natom,
589 t_atoms *atoms, rvec x[], int chainnum, gmx_bool bChange)
594 char anr[12], anm[12], anm_copy[12], altloc, resnm[12], rnr[12];
595 char xc[12], yc[12], zc[12], occup[12], bfac[12];
598 int resnr, atomnumber;
600 if (natom >= atoms->nr)
602 gmx_fatal(FARGS, "\nFound more atoms (%d) in pdb file than expected (%d)",
608 for (k = 0; (k < 5); k++, j++)
615 for (k = 0; (k < 4); k++, j++)
620 strcpy(anm_copy, anm);
626 for (k = 0; (k < 4); k++, j++)
636 for (k = 0; (k < 4); k++, j++)
642 resnr = strtol(rnr, NULL, 10);
646 /* X,Y,Z Coordinate */
647 for (k = 0; (k < 8); k++, j++)
652 for (k = 0; (k < 8); k++, j++)
657 for (k = 0; (k < 8); k++, j++)
664 for (k = 0; (k < 6); k++, j++)
671 for (k = 0; (k < 7); k++, j++)
679 atomn = &(atoms->atom[natom]);
681 atoms->resinfo[atoms->atom[natom-1].resind].nr != resnr ||
682 atoms->resinfo[atoms->atom[natom-1].resind].ic != resic ||
683 (strcmp(*atoms->resinfo[atoms->atom[natom-1].resind].name, resnm) != 0))
691 atomn->resind = atoms->atom[natom-1].resind + 1;
693 atoms->nres = atomn->resind + 1;
694 t_atoms_set_resinfo(atoms, natom, symtab, resnm, resnr, resic, chainnum, chainid);
698 atomn->resind = atoms->atom[natom-1].resind;
702 xlate_atomname_pdb2gmx(anm);
704 atoms->atomname[natom] = put_symtab(symtab, anm);
707 atomn->atomnumber = atomnumber;
708 atomn->elem[0] = '\0';
710 x[natom][XX] = strtod(xc, NULL)*0.1;
711 x[natom][YY] = strtod(yc, NULL)*0.1;
712 x[natom][ZZ] = strtod(zc, NULL)*0.1;
715 atoms->pdbinfo[natom].type = type;
716 atoms->pdbinfo[natom].atomnr = strtol(anr, NULL, 10);
717 atoms->pdbinfo[natom].altloc = altloc;
718 strcpy(atoms->pdbinfo[natom].atomnm, anm_copy);
719 atoms->pdbinfo[natom].bfac = strtod(bfac, NULL);
720 atoms->pdbinfo[natom].occup = strtod(occup, NULL);
727 gmx_bool is_hydrogen(const char *nm)
738 else if ((isdigit(buf[0])) && (buf[1] == 'H'))
745 gmx_bool is_dummymass(const char *nm)
752 if ((buf[0] == 'M') && isdigit(buf[strlen(buf)-1]))
760 static void gmx_conect_addline(gmx_conect_t *con, char *line)
763 char format[32], form2[32];
765 sprintf(form2, "%%*s");
766 sprintf(format, "%s%%d", form2);
767 if (sscanf(line, format, &ai) == 1)
771 strcat(form2, "%*s");
772 sprintf(format, "%s%%d", form2);
773 n = sscanf(line, format, &aj);
776 srenew(con->conect, ++con->nconect);
777 con->conect[con->nconect-1].ai = ai-1;
778 con->conect[con->nconect-1].aj = aj-1;
785 void gmx_conect_dump(FILE *fp, gmx_conect conect)
787 gmx_conect_t *gc = (gmx_conect_t *)conect;
790 for (i = 0; (i < gc->nconect); i++)
792 fprintf(fp, "%6s%5d%5d\n", "CONECT",
793 gc->conect[i].ai+1, gc->conect[i].aj+1);
797 gmx_conect gmx_conect_init()
803 return (gmx_conect) gc;
806 void gmx_conect_done(gmx_conect conect)
808 gmx_conect_t *gc = (gmx_conect_t *)conect;
813 gmx_bool gmx_conect_exist(gmx_conect conect, int ai, int aj)
815 gmx_conect_t *gc = (gmx_conect_t *)conect;
821 for (i = 0; (i < gc->nconect); i++)
823 if (((gc->conect[i].ai == ai) &&
824 (gc->conect[i].aj == aj)) ||
825 ((gc->conect[i].aj == ai) &&
826 (gc->conect[i].ai == aj)))
834 void gmx_conect_add(gmx_conect conect, int ai, int aj)
836 gmx_conect_t *gc = (gmx_conect_t *)conect;
842 if (!gmx_conect_exist(conect, ai, aj))
844 srenew(gc->conect, ++gc->nconect);
845 gc->conect[gc->nconect-1].ai = ai;
846 gc->conect[gc->nconect-1].aj = aj;
850 int read_pdbfile(FILE *in, char *title, int *model_nr,
851 t_atoms *atoms, rvec x[], int *ePBC, matrix box, gmx_bool bChange,
854 gmx_conect_t *gc = (gmx_conect_t *)conect;
857 gmx_bool bConnWarn = FALSE;
861 int natom, chainnum, nres_ter_prev = 0;
863 gmx_bool bStop = FALSE;
867 /* Only assume pbc when there is a CRYST1 entry */
875 open_symtab(&symtab);
881 while (!bStop && (fgets2(line, STRLEN, in) != NULL))
883 line_type = line2type(line);
889 natom = read_atom(&symtab, line, line_type, natom, atoms, x, chainnum, bChange);
895 read_anisou(line, natom, atoms);
900 read_cryst1(line, ePBC, box);
905 if (strlen(line) > 6)
908 /* skip HEADER or TITLE and spaces */
917 /* truncate after title */
931 if ((!strstr(line, ": ")) || (strstr(line+6, "MOLECULE:")))
933 if (!(c = strstr(line+6, "MOLECULE:")) )
937 /* skip 'MOLECULE:' and spaces */
946 /* truncate after title */
950 while ( (d[-1] == ';') && d > c)
979 sscanf(line, "%*s%d", model_nr);
989 gmx_conect_addline(gc, line);
993 fprintf(stderr, "WARNING: all CONECT records are ignored\n");
1003 free_symtab(&symtab);
1007 void get_pdb_coordnum(FILE *in, int *natoms)
1012 while (fgets2(line, STRLEN, in))
1014 if (strncmp(line, "ENDMDL", 6) == 0)
1018 if ((strncmp(line, "ATOM ", 6) == 0) || (strncmp(line, "HETATM", 6) == 0))
1025 void read_pdb_conf(const char *infile, char *title,
1026 t_atoms *atoms, rvec x[], int *ePBC, matrix box, gmx_bool bChange,
1031 in = gmx_fio_fopen(infile, "r");
1032 read_pdbfile(in, title, NULL, atoms, x, ePBC, box, bChange, conect);
1036 gmx_conect gmx_conect_generate(t_topology *top)
1041 /* Fill the conect records */
1042 gc = gmx_conect_init();
1044 for (f = 0; (f < F_NRE); f++)
1048 for (i = 0; (i < top->idef.il[f].nr); i += interaction_function[f].nratoms+1)
1050 gmx_conect_add(gc, top->idef.il[f].iatoms[i+1],
1051 top->idef.il[f].iatoms[i+2]);
1058 const char* get_pdbformat()
1060 static const char *pdbformat = "%-6s%5u %-4.4s%3.3s %c%4d%c %8.3f%8.3f%8.3f";
1064 const char* get_pdbformat4()
1066 static const char *pdbformat4 = "%-6s%5u %-4.4s %3.3s %c%4d%c %8.3f%8.3f%8.3f";