Version bump after 4.6.6 release
[gromacs.git] / src / gmxlib / pdbio.c
1 /*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team,
6 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012,2013, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
11 *
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
16 *
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
21 *
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 *
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
34 *
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
37 */
38 #ifdef HAVE_CONFIG_H
39 #include <config.h>
40 #endif
41
42 #include <ctype.h>
43
44 #include "sysstuff.h"
45 #include "string2.h"
46 #include "vec.h"
47 #include "smalloc.h"
48 #include "typedefs.h"
49 #include "symtab.h"
50 #include "pdbio.h"
51 #include "vec.h"
52 #include "copyrite.h"
53 #include "futil.h"
54 #include "atomprop.h"
55 #include "physics.h"
56 #include "pbc.h"
57 #include "gmxfio.h"
58
59 typedef struct {
60 int ai, aj;
61 } gmx_conection_t;
62
63 typedef struct gmx_conect_t {
64 int nconect;
65 gmx_bool bSorted;
66 gmx_conection_t *conect;
67 } gmx_conect_t;
68
69 static const char *pdbtp[epdbNR] = {
70 "ATOM ", "HETATM", "ANISOU", "CRYST1",
71 "COMPND", "MODEL", "ENDMDL", "TER", "HEADER", "TITLE", "REMARK",
72 "CONECT"
73 };
74
75
76 /* this is not very good,
77 but these are only used in gmx_trjconv and gmx_editconv */
78 static gmx_bool bWideFormat = FALSE;
79 #define REMARK_SIM_BOX "REMARK THIS IS A SIMULATION BOX"
80
81 void set_pdb_wide_format(gmx_bool bSet)
82 {
83 bWideFormat = bSet;
84 }
85
86 static void xlate_atomname_pdb2gmx(char *name)
87 {
88 int i, length;
89 char temp;
90
91 length = strlen(name);
92 if (length > 3 && isdigit(name[0]))
93 {
94 temp = name[0];
95 for (i = 1; i < length; i++)
96 {
97 name[i-1] = name[i];
98 }
99 name[length-1] = temp;
100 }
101 }
102
103 static void xlate_atomname_gmx2pdb(char *name)
104 {
105 int i, length;
106 char temp;
107
108 length = strlen(name);
109 if (length > 3 && isdigit(name[length-1]))
110 {
111 temp = name[length-1];
112 for (i = length-1; i > 0; --i)
113 {
114 name[i] = name[i-1];
115 }
116 name[0] = temp;
117 }
118 }
119
120
121 void gmx_write_pdb_box(FILE *out, int ePBC, matrix box)
122 {
123 real alpha, beta, gamma;
124
125 if (ePBC == -1)
126 {
127 ePBC = guess_ePBC(box);
128 }
129
130 if (ePBC == epbcNONE)
131 {
132 return;
133 }
134
135 if (norm2(box[YY])*norm2(box[ZZ]) != 0)
136 {
137 alpha = RAD2DEG*acos(cos_angle_no_table(box[YY], box[ZZ]));
138 }
139 else
140 {
141 alpha = 90;
142 }
143 if (norm2(box[XX])*norm2(box[ZZ]) != 0)
144 {
145 beta = RAD2DEG*acos(cos_angle_no_table(box[XX], box[ZZ]));
146 }
147 else
148 {
149 beta = 90;
150 }
151 if (norm2(box[XX])*norm2(box[YY]) != 0)
152 {
153 gamma = RAD2DEG*acos(cos_angle_no_table(box[XX], box[YY]));
154 }
155 else
156 {
157 gamma = 90;
158 }
159 fprintf(out, "REMARK THIS IS A SIMULATION BOX\n");
160 if (ePBC != epbcSCREW)
161 {
162 fprintf(out, "CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f %-11s%4d\n",
163 10*norm(box[XX]), 10*norm(box[YY]), 10*norm(box[ZZ]),
164 alpha, beta, gamma, "P 1", 1);
165 }
166 else
167 {
168 /* Double the a-vector length and write the correct space group */
169 fprintf(out, "CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f %-11s%4d\n",
170 20*norm(box[XX]), 10*norm(box[YY]), 10*norm(box[ZZ]),
171 alpha, beta, gamma, "P 21 1 1", 1);
172
173 }
174 }
175
176 static void read_cryst1(char *line, int *ePBC, matrix box)
177 {
178 #define SG_SIZE 11
179 char sa[12], sb[12], sc[12], sg[SG_SIZE+1], ident;
180 double fa, fb, fc, alpha, beta, gamma, cosa, cosb, cosg, sing;
181 int syma, symb, symc;
182 int ePBC_file;
183
184 sscanf(line, "%*s%s%s%s%lf%lf%lf", sa, sb, sc, &alpha, &beta, &gamma);
185
186 ePBC_file = -1;
187 if (strlen(line) >= 55)
188 {
189 strncpy(sg, line+55, SG_SIZE);
190 sg[SG_SIZE] = '\0';
191 ident = ' ';
192 syma = 0;
193 symb = 0;
194 symc = 0;
195 sscanf(sg, "%c %d %d %d", &ident, &syma, &symb, &symc);
196 if (ident == 'P' && syma == 1 && symb <= 1 && symc <= 1)
197 {
198 fc = strtod(sc, NULL)*0.1;
199 ePBC_file = (fc > 0 ? epbcXYZ : epbcXY);
200 }
201 if (ident == 'P' && syma == 21 && symb == 1 && symc == 1)
202 {
203 ePBC_file = epbcSCREW;
204 }
205 }
206 if (ePBC)
207 {
208 *ePBC = ePBC_file;
209 }
210
211 if (box)
212 {
213 fa = strtod(sa, NULL)*0.1;
214 fb = strtod(sb, NULL)*0.1;
215 fc = strtod(sc, NULL)*0.1;
216 if (ePBC_file == epbcSCREW)
217 {
218 fa *= 0.5;
219 }
220 clear_mat(box);
221 box[XX][XX] = fa;
222 if ((alpha != 90.0) || (beta != 90.0) || (gamma != 90.0))
223 {
224 if (alpha != 90.0)
225 {
226 cosa = cos(alpha*DEG2RAD);
227 }
228 else
229 {
230 cosa = 0;
231 }
232 if (beta != 90.0)
233 {
234 cosb = cos(beta*DEG2RAD);
235 }
236 else
237 {
238 cosb = 0;
239 }
240 if (gamma != 90.0)
241 {
242 cosg = cos(gamma*DEG2RAD);
243 sing = sin(gamma*DEG2RAD);
244 }
245 else
246 {
247 cosg = 0;
248 sing = 1;
249 }
250 box[YY][XX] = fb*cosg;
251 box[YY][YY] = fb*sing;
252 box[ZZ][XX] = fc*cosb;
253 box[ZZ][YY] = fc*(cosa - cosb*cosg)/sing;
254 box[ZZ][ZZ] = sqrt(fc*fc
255 - box[ZZ][XX]*box[ZZ][XX] - box[ZZ][YY]*box[ZZ][YY]);
256 }
257 else
258 {
259 box[YY][YY] = fb;
260 box[ZZ][ZZ] = fc;
261 }
262 }
263 }
264
265 void write_pdbfile_indexed(FILE *out, const char *title,
266 t_atoms *atoms, rvec x[],
267 int ePBC, matrix box, char chainid,
268 int model_nr, atom_id nindex, atom_id index[],
269 gmx_conect conect, gmx_bool bTerSepChains)
270 {
271 gmx_conect_t *gc = (gmx_conect_t *)conect;
272 char resnm[6], nm[6], pdbform[128], pukestring[100];
273 atom_id i, ii;
274 int resind, resnr, type;
275 unsigned char resic, ch;
276 real occup, bfac;
277 gmx_bool bOccup;
278 int nlongname = 0;
279 int chainnum, lastchainnum;
280 int lastresind, lastchainresind;
281 gmx_residuetype_t rt;
282 const char *p_restype;
283 const char *p_lastrestype;
284
285 gmx_residuetype_init(&rt);
286
287 bromacs(pukestring, 99);
288 fprintf(out, "TITLE %s\n", (title && title[0]) ? title : pukestring);
289 if (bWideFormat)
290 {
291 fprintf(out, "REMARK This file does not adhere to the PDB standard\n");
292 fprintf(out, "REMARK As a result of, some programs may not like it\n");
293 }
294 if (box && ( norm2(box[XX]) || norm2(box[YY]) || norm2(box[ZZ]) ) )
295 {
296 gmx_write_pdb_box(out, ePBC, box);
297 }
298 if (atoms->pdbinfo)
299 {
300 /* Check whether any occupancies are set, in that case leave it as is,
301 * otherwise set them all to one
302 */
303 bOccup = TRUE;
304 for (ii = 0; (ii < nindex) && bOccup; ii++)
305 {
306 i = index[ii];
307 bOccup = bOccup && (atoms->pdbinfo[i].occup == 0.0);
308 }
309 }
310 else
311 {
312 bOccup = FALSE;
313 }
314
315 fprintf(out, "MODEL %8d\n", model_nr > 0 ? model_nr : 1);
316
317 lastchainresind = -1;
318 lastchainnum = -1;
319 resind = -1;
320 p_restype = NULL;
321
322 for (ii = 0; ii < nindex; ii++)
323 {
324 i = index[ii];
325 lastresind = resind;
326 resind = atoms->atom[i].resind;
327 chainnum = atoms->resinfo[resind].chainnum;
328 p_lastrestype = p_restype;
329 gmx_residuetype_get_type(rt, *atoms->resinfo[resind].name, &p_restype);
330
331 /* Add a TER record if we changed chain, and if either the previous or this chain is protein/DNA/RNA. */
332 if (bTerSepChains && ii > 0 && chainnum != lastchainnum)
333 {
334 /* Only add TER if the previous chain contained protein/DNA/RNA. */
335 if (gmx_residuetype_is_protein(rt, p_lastrestype) || gmx_residuetype_is_dna(rt, p_lastrestype) || gmx_residuetype_is_rna(rt, p_lastrestype))
336 {
337 fprintf(out, "TER\n");
338 }
339 lastchainnum = chainnum;
340 lastchainresind = lastresind;
341 }
342
343 strncpy(resnm, *atoms->resinfo[resind].name, sizeof(resnm)-1);
344 strncpy(nm, *atoms->atomname[i], sizeof(nm)-1);
345 /* rename HG12 to 2HG1, etc. */
346 xlate_atomname_gmx2pdb(nm);
347 resnr = atoms->resinfo[resind].nr;
348 resic = atoms->resinfo[resind].ic;
349 if (chainid != ' ')
350 {
351 ch = chainid;
352 }
353 else
354 {
355 ch = atoms->resinfo[resind].chainid;
356
357 if (ch == 0)
358 {
359 ch = ' ';
360 }
361 }
362 if (resnr >= 10000)
363 {
364 resnr = resnr % 10000;
365 }
366 if (atoms->pdbinfo)
367 {
368 type = atoms->pdbinfo[i].type;
369 occup = bOccup ? 1.0 : atoms->pdbinfo[i].occup;
370 bfac = atoms->pdbinfo[i].bfac;
371 }
372 else
373 {
374 type = 0;
375 occup = 1.0;
376 bfac = 0.0;
377 }
378 if (bWideFormat)
379 {
380 strcpy(pdbform,
381 "%-6s%5u %-4.4s %3.3s %c%4d%c %10.5f%10.5f%10.5f%8.4f%8.4f %2s\n");
382 }
383 else
384 {
385 /* Check whether atomname is an element name */
386 if ((strlen(nm) < 4) && (gmx_strcasecmp(nm, atoms->atom[i].elem) != 0))
387 {
388 strcpy(pdbform, pdbformat);
389 }
390 else
391 {
392 strcpy(pdbform, pdbformat4);
393 if (strlen(nm) > 4)
394 {
395 int maxwln = 20;
396 if (nlongname < maxwln)
397 {
398 fprintf(stderr, "WARNING: Writing out atom name (%s) longer than 4 characters to .pdb file\n", nm);
399 }
400 else if (nlongname == maxwln)
401 {
402 fprintf(stderr, "WARNING: More than %d long atom names, will not write more warnings\n", maxwln);
403 }
404 nlongname++;
405 }
406 }
407 strcat(pdbform, "%6.2f%6.2f %2s\n");
408 }
409 fprintf(out, pdbform, pdbtp[type], (i+1)%100000, nm, resnm, ch, resnr,
410 (resic == '\0') ? ' ' : resic,
411 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], occup, bfac, atoms->atom[i].elem);
412 if (atoms->pdbinfo && atoms->pdbinfo[i].bAnisotropic)
413 {
414 fprintf(out, "ANISOU%5u %-4.4s%3.3s %c%4d%c %7d%7d%7d%7d%7d%7d\n",
415 (i+1)%100000, nm, resnm, ch, resnr,
416 (resic == '\0') ? ' ' : resic,
417 atoms->pdbinfo[i].uij[0], atoms->pdbinfo[i].uij[1],
418 atoms->pdbinfo[i].uij[2], atoms->pdbinfo[i].uij[3],
419 atoms->pdbinfo[i].uij[4], atoms->pdbinfo[i].uij[5]);
420 }
421 }
422
423 fprintf(out, "TER\n");
424 fprintf(out, "ENDMDL\n");
425
426 if (NULL != gc)
427 {
428 /* Write conect records */
429 for (i = 0; (i < gc->nconect); i++)
430 {
431 fprintf(out, "CONECT%5d%5d\n", gc->conect[i].ai+1, gc->conect[i].aj+1);
432 }
433 }
434
435 gmx_residuetype_destroy(rt);
436 }
437
438 void write_pdbfile(FILE *out, const char *title, t_atoms *atoms, rvec x[],
439 int ePBC, matrix box, char chainid, int model_nr, gmx_conect conect, gmx_bool bTerSepChains)
440 {
441 atom_id i, *index;
442
443 snew(index, atoms->nr);
444 for (i = 0; i < atoms->nr; i++)
445 {
446 index[i] = i;
447 }
448 write_pdbfile_indexed(out, title, atoms, x, ePBC, box, chainid, model_nr,
449 atoms->nr, index, conect, bTerSepChains);
450 sfree(index);
451 }
452
453 int line2type(char *line)
454 {
455 int k;
456 char type[8];
457
458 for (k = 0; (k < 6); k++)
459 {
460 type[k] = line[k];
461 }
462 type[k] = '\0';
463
464 for (k = 0; (k < epdbNR); k++)
465 {
466 if (strncmp(type, pdbtp[k], strlen(pdbtp[k])) == 0)
467 {
468 break;
469 }
470 }
471
472 return k;
473 }
474
475 static void read_anisou(char line[], int natom, t_atoms *atoms)
476 {
477 int i, j, k, atomnr;
478 char nc = '\0';
479 char anr[12], anm[12];
480
481 /* Skip over type */
482 j = 6;
483 for (k = 0; (k < 5); k++, j++)
484 {
485 anr[k] = line[j];
486 }
487 anr[k] = nc;
488 j++;
489 for (k = 0; (k < 4); k++, j++)
490 {
491 anm[k] = line[j];
492 }
493 anm[k] = nc;
494 j++;
495
496 /* Strip off spaces */
497 trim(anm);
498
499 /* Search backwards for number and name only */
500 atomnr = strtol(anr, NULL, 10);
501 for (i = natom-1; (i >= 0); i--)
502 {
503 if ((strcmp(anm, *(atoms->atomname[i])) == 0) &&
504 (atomnr == atoms->pdbinfo[i].atomnr))
505 {
506 break;
507 }
508 }
509 if (i < 0)
510 {
511 fprintf(stderr, "Skipping ANISOU record (atom %s %d not found)\n",
512 anm, atomnr);
513 }
514 else
515 {
516 if (sscanf(line+29, "%d%d%d%d%d%d",
517 &atoms->pdbinfo[i].uij[U11], &atoms->pdbinfo[i].uij[U22],
518 &atoms->pdbinfo[i].uij[U33], &atoms->pdbinfo[i].uij[U12],
519 &atoms->pdbinfo[i].uij[U13], &atoms->pdbinfo[i].uij[U23])
520 == 6)
521 {
522 atoms->pdbinfo[i].bAnisotropic = TRUE;
523 }
524 else
525 {
526 fprintf(stderr, "Invalid ANISOU record for atom %d\n", i);
527 atoms->pdbinfo[i].bAnisotropic = FALSE;
528 }
529 }
530 }
531
532 void get_pdb_atomnumber(t_atoms *atoms, gmx_atomprop_t aps)
533 {
534 int i, atomnumber, len;
535 size_t k;
536 char anm[6], anm_copy[6], *ptr;
537 char nc = '\0';
538 real eval;
539
540 if (!atoms->pdbinfo)
541 {
542 gmx_incons("Trying to deduce atomnumbers when no pdb information is present");
543 }
544 for (i = 0; (i < atoms->nr); i++)
545 {
546 strcpy(anm, atoms->pdbinfo[i].atomnm);
547 strcpy(anm_copy, atoms->pdbinfo[i].atomnm);
548 len = strlen(anm);
549 atomnumber = NOTSET;
550 if ((anm[0] != ' ') && ((len <= 2) || ((len > 2) && !isdigit(anm[2]))))
551 {
552 anm_copy[2] = nc;
553 if (gmx_atomprop_query(aps, epropElement, "???", anm_copy, &eval))
554 {
555 atomnumber = gmx_nint(eval);
556 }
557 else
558 {
559 anm_copy[1] = nc;
560 if (gmx_atomprop_query(aps, epropElement, "???", anm_copy, &eval))
561 {
562 atomnumber = gmx_nint(eval);
563 }
564 }
565 }
566 if (atomnumber == NOTSET)
567 {
568 k = 0;
569 while ((k < strlen(anm)) && (isspace(anm[k]) || isdigit(anm[k])))
570 {
571 k++;
572 }
573 anm_copy[0] = anm[k];
574 anm_copy[1] = nc;
575 if (gmx_atomprop_query(aps, epropElement, "???", anm_copy, &eval))
576 {
577 atomnumber = gmx_nint(eval);
578 }
579 }
580 atoms->atom[i].atomnumber = atomnumber;
581 ptr = gmx_atomprop_element(aps, atomnumber);
582 strncpy(atoms->atom[i].elem, ptr == NULL ? "" : ptr, 4);
583 if (debug)
584 {
585 fprintf(debug, "Atomnumber for atom '%s' is %d\n", anm, atomnumber);
586 }
587 }
588 }
589
590 static int read_atom(t_symtab *symtab,
591 char line[], int type, int natom,
592 t_atoms *atoms, rvec x[], int chainnum, gmx_bool bChange)
593 {
594 t_atom *atomn;
595 int j, k;
596 char nc = '\0';
597 char anr[12], anm[12], anm_copy[12], altloc, resnm[12], rnr[12];
598 char xc[12], yc[12], zc[12], occup[12], bfac[12];
599 unsigned char resic;
600 char chainid;
601 int resnr, atomnumber;
602
603 if (natom >= atoms->nr)
604 {
605 gmx_fatal(FARGS, "\nFound more atoms (%d) in pdb file than expected (%d)",
606 natom+1, atoms->nr);
607 }
608
609 /* Skip over type */
610 j = 6;
611 for (k = 0; (k < 5); k++, j++)
612 {
613 anr[k] = line[j];
614 }
615 anr[k] = nc;
616 trim(anr);
617 j++;
618 for (k = 0; (k < 4); k++, j++)
619 {
620 anm[k] = line[j];
621 }
622 anm[k] = nc;
623 strcpy(anm_copy, anm);
624 rtrim(anm_copy);
625 atomnumber = NOTSET;
626 trim(anm);
627 altloc = line[j];
628 j++;
629 for (k = 0; (k < 4); k++, j++)
630 {
631 resnm[k] = line[j];
632 }
633 resnm[k] = nc;
634 trim(resnm);
635
636 chainid = line[j];
637 j++;
638
639 for (k = 0; (k < 4); k++, j++)
640 {
641 rnr[k] = line[j];
642 }
643 rnr[k] = nc;
644 trim(rnr);
645 resnr = strtol(rnr, NULL, 10);
646 resic = line[j];
647 j += 4;
648
649 /* X,Y,Z Coordinate */
650 for (k = 0; (k < 8); k++, j++)
651 {
652 xc[k] = line[j];
653 }
654 xc[k] = nc;
655 for (k = 0; (k < 8); k++, j++)
656 {
657 yc[k] = line[j];
658 }
659 yc[k] = nc;
660 for (k = 0; (k < 8); k++, j++)
661 {
662 zc[k] = line[j];
663 }
664 zc[k] = nc;
665
666 /* Weight */
667 for (k = 0; (k < 6); k++, j++)
668 {
669 occup[k] = line[j];
670 }
671 occup[k] = nc;
672
673 /* B-Factor */
674 for (k = 0; (k < 7); k++, j++)
675 {
676 bfac[k] = line[j];
677 }
678 bfac[k] = nc;
679
680 if (atoms->atom)
681 {
682 atomn = &(atoms->atom[natom]);
683 if ((natom == 0) ||
684 atoms->resinfo[atoms->atom[natom-1].resind].nr != resnr ||
685 atoms->resinfo[atoms->atom[natom-1].resind].ic != resic ||
686 (strcmp(*atoms->resinfo[atoms->atom[natom-1].resind].name, resnm) != 0))
687 {
688 if (natom == 0)
689 {
690 atomn->resind = 0;
691 }
692 else
693 {
694 atomn->resind = atoms->atom[natom-1].resind + 1;
695 }
696 atoms->nres = atomn->resind + 1;
697 t_atoms_set_resinfo(atoms, natom, symtab, resnm, resnr, resic, chainnum, chainid);
698 }
699 else
700 {
701 atomn->resind = atoms->atom[natom-1].resind;
702 }
703 if (bChange)
704 {
705 xlate_atomname_pdb2gmx(anm);
706 }
707 atoms->atomname[natom] = put_symtab(symtab, anm);
708 atomn->m = 0.0;
709 atomn->q = 0.0;
710 atomn->atomnumber = atomnumber;
711 atomn->elem[0] = '\0';
712 }
713 x[natom][XX] = strtod(xc, NULL)*0.1;
714 x[natom][YY] = strtod(yc, NULL)*0.1;
715 x[natom][ZZ] = strtod(zc, NULL)*0.1;
716 if (atoms->pdbinfo)
717 {
718 atoms->pdbinfo[natom].type = type;
719 atoms->pdbinfo[natom].atomnr = strtol(anr, NULL, 10);
720 atoms->pdbinfo[natom].altloc = altloc;
721 strcpy(atoms->pdbinfo[natom].atomnm, anm_copy);
722 atoms->pdbinfo[natom].bfac = strtod(bfac, NULL);
723 atoms->pdbinfo[natom].occup = strtod(occup, NULL);
724 }
725 natom++;
726
727 return natom;
728 }
729
730 gmx_bool is_hydrogen(const char *nm)
731 {
732 char buf[30];
733
734 strcpy(buf, nm);
735 trim(buf);
736
737 if (buf[0] == 'H')
738 {
739 return TRUE;
740 }
741 else if ((isdigit(buf[0])) && (buf[1] == 'H'))
742 {
743 return TRUE;
744 }
745 return FALSE;
746 }
747
748 gmx_bool is_dummymass(const char *nm)
749 {
750 char buf[30];
751
752 strcpy(buf, nm);
753 trim(buf);
754
755 if ((buf[0] == 'M') && isdigit(buf[strlen(buf)-1]))
756 {
757 return TRUE;
758 }
759
760 return FALSE;
761 }
762
763 static void gmx_conect_addline(gmx_conect_t *con, char *line)
764 {
765 int n, ai, aj;
766 char format[32], form2[32];
767
768 sprintf(form2, "%%*s");
769 sprintf(format, "%s%%d", form2);
770 if (sscanf(line, format, &ai) == 1)
771 {
772 do
773 {
774 strcat(form2, "%*s");
775 sprintf(format, "%s%%d", form2);
776 n = sscanf(line, format, &aj);
777 if (n == 1)
778 {
779 srenew(con->conect, ++con->nconect);
780 con->conect[con->nconect-1].ai = ai-1;
781 con->conect[con->nconect-1].aj = aj-1;
782 }
783 }
784 while (n == 1);
785 }
786 }
787
788 void gmx_conect_dump(FILE *fp, gmx_conect conect)
789 {
790 gmx_conect_t *gc = (gmx_conect_t *)conect;
791 int i;
792
793 for (i = 0; (i < gc->nconect); i++)
794 {
795 fprintf(fp, "%6s%5d%5d\n", "CONECT",
796 gc->conect[i].ai+1, gc->conect[i].aj+1);
797 }
798 }
799
800 gmx_conect gmx_conect_init()
801 {
802 gmx_conect_t *gc;
803
804 snew(gc, 1);
805
806 return (gmx_conect) gc;
807 }
808
809 void gmx_conect_done(gmx_conect conect)
810 {
811 gmx_conect_t *gc = (gmx_conect_t *)conect;
812
813 sfree(gc->conect);
814 }
815
816 gmx_bool gmx_conect_exist(gmx_conect conect, int ai, int aj)
817 {
818 gmx_conect_t *gc = (gmx_conect_t *)conect;
819 int i;
820
821 /* if (!gc->bSorted)
822 sort_conect(gc);*/
823
824 for (i = 0; (i < gc->nconect); i++)
825 {
826 if (((gc->conect[i].ai == ai) &&
827 (gc->conect[i].aj == aj)) ||
828 ((gc->conect[i].aj == ai) &&
829 (gc->conect[i].ai == aj)))
830 {
831 return TRUE;
832 }
833 }
834 return FALSE;
835 }
836
837 void gmx_conect_add(gmx_conect conect, int ai, int aj)
838 {
839 gmx_conect_t *gc = (gmx_conect_t *)conect;
840 int i;
841
842 /* if (!gc->bSorted)
843 sort_conect(gc);*/
844
845 if (!gmx_conect_exist(conect, ai, aj))
846 {
847 srenew(gc->conect, ++gc->nconect);
848 gc->conect[gc->nconect-1].ai = ai;
849 gc->conect[gc->nconect-1].aj = aj;
850 }
851 }
852
853 int read_pdbfile(FILE *in, char *title, int *model_nr,
854 t_atoms *atoms, rvec x[], int *ePBC, matrix box, gmx_bool bChange,
855 gmx_conect conect)
856 {
857 gmx_conect_t *gc = (gmx_conect_t *)conect;
858 t_symtab symtab;
859 gmx_bool bCOMPND;
860 gmx_bool bConnWarn = FALSE;
861 char line[STRLEN+1];
862 int line_type;
863 char *c, *d;
864 int natom, chainnum, nres_ter_prev = 0;
865 char chidmax = ' ';
866 gmx_bool bStop = FALSE;
867
868 if (ePBC)
869 {
870 /* Only assume pbc when there is a CRYST1 entry */
871 *ePBC = epbcNONE;
872 }
873 if (box != NULL)
874 {
875 clear_mat(box);
876 }
877
878 open_symtab(&symtab);
879
880 bCOMPND = FALSE;
881 title[0] = '\0';
882 natom = 0;
883 chainnum = 0;
884 while (!bStop && (fgets2(line, STRLEN, in) != NULL))
885 {
886 line_type = line2type(line);
887
888 switch (line_type)
889 {
890 case epdbATOM:
891 case epdbHETATM:
892 natom = read_atom(&symtab, line, line_type, natom, atoms, x, chainnum, bChange);
893 break;
894
895 case epdbANISOU:
896 if (atoms->pdbinfo)
897 {
898 read_anisou(line, natom, atoms);
899 }
900 break;
901
902 case epdbCRYST1:
903 read_cryst1(line, ePBC, box);
904 break;
905
906 case epdbTITLE:
907 case epdbHEADER:
908 if (strlen(line) > 6)
909 {
910 c = line+6;
911 /* skip HEADER or TITLE and spaces */
912 while (c[0] != ' ')
913 {
914 c++;
915 }
916 while (c[0] == ' ')
917 {
918 c++;
919 }
920 /* truncate after title */
921 d = strstr(c, " ");
922 if (d)
923 {
924 d[0] = '\0';
925 }
926 if (strlen(c) > 0)
927 {
928 strcpy(title, c);
929 }
930 }
931 break;
932
933 case epdbCOMPND:
934 if ((!strstr(line, ": ")) || (strstr(line+6, "MOLECULE:")))
935 {
936 if (!(c = strstr(line+6, "MOLECULE:")) )
937 {
938 c = line;
939 }
940 /* skip 'MOLECULE:' and spaces */
941 while (c[0] != ' ')
942 {
943 c++;
944 }
945 while (c[0] == ' ')
946 {
947 c++;
948 }
949 /* truncate after title */
950 d = strstr(c, " ");
951 if (d)
952 {
953 while ( (d[-1] == ';') && d > c)
954 {
955 d--;
956 }
957 d[0] = '\0';
958 }
959 if (strlen(c) > 0)
960 {
961 if (bCOMPND)
962 {
963 strcat(title, "; ");
964 strcat(title, c);
965 }
966 else
967 {
968 strcpy(title, c);
969 }
970 }
971 bCOMPND = TRUE;
972 }
973 break;
974
975 case epdbTER:
976 chainnum++;
977 break;
978
979 case epdbMODEL:
980 if (model_nr)
981 {
982 sscanf(line, "%*s%d", model_nr);
983 }
984 break;
985
986 case epdbENDMDL:
987 bStop = TRUE;
988 break;
989 case epdbCONECT:
990 if (gc)
991 {
992 gmx_conect_addline(gc, line);
993 }
994 else if (!bConnWarn)
995 {
996 fprintf(stderr, "WARNING: all CONECT records are ignored\n");
997 bConnWarn = TRUE;
998 }
999 break;
1000
1001 default:
1002 break;
1003 }
1004 }
1005
1006 free_symtab(&symtab);
1007 return natom;
1008 }
1009
1010 void get_pdb_coordnum(FILE *in, int *natoms)
1011 {
1012 char line[STRLEN];
1013
1014 *natoms = 0;
1015 while (fgets2(line, STRLEN, in))
1016 {
1017 if (strncmp(line, "ENDMDL", 6) == 0)
1018 {
1019 break;
1020 }
1021 if ((strncmp(line, "ATOM ", 6) == 0) || (strncmp(line, "HETATM", 6) == 0))
1022 {
1023 (*natoms)++;
1024 }
1025 }
1026 }
1027
1028 void read_pdb_conf(const char *infile, char *title,
1029 t_atoms *atoms, rvec x[], int *ePBC, matrix box, gmx_bool bChange,
1030 gmx_conect conect)
1031 {
1032 FILE *in;
1033
1034 in = gmx_fio_fopen(infile, "r");
1035 read_pdbfile(in, title, NULL, atoms, x, ePBC, box, bChange, conect);
1036 gmx_fio_fclose(in);
1037 }
1038
1039 gmx_conect gmx_conect_generate(t_topology *top)
1040 {
1041 int f, i;
1042 gmx_conect gc;
1043
1044 /* Fill the conect records */
1045 gc = gmx_conect_init();
1046
1047 for (f = 0; (f < F_NRE); f++)
1048 {
1049 if (IS_CHEMBOND(f))
1050 {
1051 for (i = 0; (i < top->idef.il[f].nr); i += interaction_function[f].nratoms+1)
1052 {
1053 gmx_conect_add(gc, top->idef.il[f].iatoms[i+1],
1054 top->idef.il[f].iatoms[i+2]);
1055 }
1056 }
1057 }
1058 return gc;
1059 }