?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Add readConfAndTopology
[gromacs.git] / src / gromacs / fileio / pdbio.cpp
blob7824afd521cc57389ecb8c2a5d23b80b389ed313
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 * Copyright (c) 2013,2014,2015,2016, 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.
10 *
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.
15 *
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.
20 *
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.
25 *
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.
33 *
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.
36 */
37 #include "gmxpre.h"
38
39 #include "pdbio.h"
40
41 #include <cctype>
42 #include <cmath>
43 #include <cstdlib>
44 #include <cstring>
45
46 #include <string>
47
48 #include "gromacs/fileio/gmxfio.h"
49 #include "gromacs/math/units.h"
50 #include "gromacs/math/vec.h"
51 #include "gromacs/pbcutil/pbc.h"
52 #include "gromacs/topology/atomprop.h"
53 #include "gromacs/topology/ifunc.h"
54 #include "gromacs/topology/residuetypes.h"
55 #include "gromacs/topology/symtab.h"
56 #include "gromacs/topology/topology.h"
57 #include "gromacs/utility/coolstuff.h"
58 #include "gromacs/utility/cstringutil.h"
59 #include "gromacs/utility/fatalerror.h"
60 #include "gromacs/utility/futil.h"
61 #include "gromacs/utility/smalloc.h"
62 #include "gromacs/utility/snprintf.h"
63
64 typedef struct {
65 int ai, aj;
66 } gmx_conection_t;
67
68 typedef struct gmx_conect_t {
69 int nconect;
70 gmx_bool bSorted;
71 gmx_conection_t *conect;
72 } gmx_conect_t;
73
74 static const char *pdbtp[epdbNR] = {
75 "ATOM ", "HETATM", "ANISOU", "CRYST1",
76 "COMPND", "MODEL", "ENDMDL", "TER", "HEADER", "TITLE", "REMARK",
77 "CONECT"
78 };
79
80
81 #define REMARK_SIM_BOX "REMARK THIS IS A SIMULATION BOX"
82
83 static void xlate_atomname_pdb2gmx(char *name)
84 {
85 int i, length;
86 char temp;
87
88 length = std::strlen(name);
89 if (length > 3 && std::isdigit(name[0]))
90 {
91 temp = name[0];
92 for (i = 1; i < length; i++)
93 {
94 name[i-1] = name[i];
95 }
96 name[length-1] = temp;
97 }
98 }
99
100 static void xlate_atomname_gmx2pdb(char *name)
101 {
102 int i, length;
103 char temp;
104
105 length = std::strlen(name);
106 if (length > 3 && std::isdigit(name[length-1]))
107 {
108 temp = name[length-1];
109 for (i = length-1; i > 0; --i)
110 {
111 name[i] = name[i-1];
112 }
113 name[0] = temp;
114 }
115 }
116
117
118 void gmx_write_pdb_box(FILE *out, int ePBC, const matrix box)
119 {
120 real alpha, beta, gamma;
121
122 if (ePBC == -1)
123 {
124 ePBC = guess_ePBC(box);
125 }
126
127 if (ePBC == epbcNONE)
128 {
129 return;
130 }
131
132 if (norm2(box[YY])*norm2(box[ZZ]) != 0)
133 {
134 alpha = RAD2DEG*gmx_angle(box[YY], box[ZZ]);
135 }
136 else
137 {
138 alpha = 90;
139 }
140 if (norm2(box[XX])*norm2(box[ZZ]) != 0)
141 {
142 beta = RAD2DEG*gmx_angle(box[XX], box[ZZ]);
143 }
144 else
145 {
146 beta = 90;
147 }
148 if (norm2(box[XX])*norm2(box[YY]) != 0)
149 {
150 gamma = RAD2DEG*gmx_angle(box[XX], box[YY]);
151 }
152 else
153 {
154 gamma = 90;
155 }
156 fprintf(out, "REMARK THIS IS A SIMULATION BOX\n");
157 if (ePBC != epbcSCREW)
158 {
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);
162 }
163 else
164 {
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);
169
170 }
171 }
172
173 static void read_cryst1(char *line, int *ePBC, matrix box)
174 {
175 #define SG_SIZE 11
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;
179 int ePBC_file;
180
181 sscanf(line, "%*s%s%s%s%lf%lf%lf", sa, sb, sc, &alpha, &beta, &gamma);
182
183 ePBC_file = -1;
184 if (strlen(line) >= 55)
185 {
186 strncpy(sg, line+55, SG_SIZE);
187 sg[SG_SIZE] = '\0';
188 ident = ' ';
189 syma = 0;
190 symb = 0;
191 symc = 0;
192 sscanf(sg, "%c %d %d %d", &ident, &syma, &symb, &symc);
193 if (ident == 'P' && syma == 1 && symb <= 1 && symc <= 1)
194 {
195 fc = strtod(sc, NULL)*0.1;
196 ePBC_file = (fc > 0 ? epbcXYZ : epbcXY);
197 }
198 if (ident == 'P' && syma == 21 && symb == 1 && symc == 1)
199 {
200 ePBC_file = epbcSCREW;
201 }
202 }
203 if (ePBC)
204 {
205 *ePBC = ePBC_file;
206 }
207
208 if (box)
209 {
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)
214 {
215 fa *= 0.5;
216 }
217 clear_mat(box);
218 box[XX][XX] = fa;
219 if ((alpha != 90.0) || (beta != 90.0) || (gamma != 90.0))
220 {
221 if (alpha != 90.0)
222 {
223 cosa = std::cos(alpha*DEG2RAD);
224 }
225 else
226 {
227 cosa = 0;
228 }
229 if (beta != 90.0)
230 {
231 cosb = std::cos(beta*DEG2RAD);
232 }
233 else
234 {
235 cosb = 0;
236 }
237 if (gamma != 90.0)
238 {
239 cosg = std::cos(gamma*DEG2RAD);
240 sing = std::sin(gamma*DEG2RAD);
241 }
242 else
243 {
244 cosg = 0;
245 sing = 1;
246 }
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] = std::sqrt(fc*fc
252 - box[ZZ][XX]*box[ZZ][XX] - box[ZZ][YY]*box[ZZ][YY]);
253 }
254 else
255 {
256 box[YY][YY] = fb;
257 box[ZZ][ZZ] = fc;
258 }
259 }
260 }
261
262 void write_pdbfile_indexed(FILE *out, const char *title,
263 const t_atoms *atoms, const rvec x[],
264 int ePBC, const matrix box, char chainid,
265 int model_nr, int nindex, const int index[],
266 gmx_conect conect, gmx_bool bTerSepChains)
267 {
268 gmx_conect_t *gc = (gmx_conect_t *)conect;
269 char resnm[6], nm[6];
270 int i, ii;
271 int resind, resnr;
272 enum PDB_record type;
273 unsigned char resic, ch;
274 char altloc;
275 real occup, bfac;
276 gmx_bool bOccup;
277 int chainnum, lastchainnum;
278 gmx_residuetype_t*rt;
279 const char *p_restype;
280 const char *p_lastrestype;
281
282 gmx_residuetype_init(&rt);
283
284 fprintf(out, "TITLE %s\n", (title && title[0]) ? title : gmx::bromacs().c_str());
285 if (box && ( norm2(box[XX]) || norm2(box[YY]) || norm2(box[ZZ]) ) )
286 {
287 gmx_write_pdb_box(out, ePBC, box);
288 }
289 if (atoms->havePdbInfo)
290 {
291 /* Check whether any occupancies are set, in that case leave it as is,
292 * otherwise set them all to one
293 */
294 bOccup = TRUE;
295 for (ii = 0; (ii < nindex) && bOccup; ii++)
296 {
297 i = index[ii];
298 bOccup = bOccup && (atoms->pdbinfo[i].occup == 0.0);
299 }
300 }
301 else
302 {
303 bOccup = FALSE;
304 }
305
306 fprintf(out, "MODEL %8d\n", model_nr > 0 ? model_nr : 1);
307
308 lastchainnum = -1;
309 p_restype = NULL;
310
311 for (ii = 0; ii < nindex; ii++)
312 {
313 i = index[ii];
314 resind = atoms->atom[i].resind;
315 chainnum = atoms->resinfo[resind].chainnum;
316 p_lastrestype = p_restype;
317 gmx_residuetype_get_type(rt, *atoms->resinfo[resind].name, &p_restype);
318
319 /* Add a TER record if we changed chain, and if either the previous or this chain is protein/DNA/RNA. */
320 if (bTerSepChains && ii > 0 && chainnum != lastchainnum)
321 {
322 /* Only add TER if the previous chain contained protein/DNA/RNA. */
323 if (gmx_residuetype_is_protein(rt, p_lastrestype) || gmx_residuetype_is_dna(rt, p_lastrestype) || gmx_residuetype_is_rna(rt, p_lastrestype))
324 {
325 fprintf(out, "TER\n");
326 }
327 lastchainnum = chainnum;
328 }
329
330 strncpy(resnm, *atoms->resinfo[resind].name, sizeof(resnm)-1);
331 resnm[sizeof(resnm)-1] = 0;
332 strncpy(nm, *atoms->atomname[i], sizeof(nm)-1);
333 nm[sizeof(nm)-1] = 0;
334
335 /* rename HG12 to 2HG1, etc. */
336 xlate_atomname_gmx2pdb(nm);
337 resnr = atoms->resinfo[resind].nr;
338 resic = atoms->resinfo[resind].ic;
339 if (chainid != ' ')
340 {
341 ch = chainid;
342 }
343 else
344 {
345 ch = atoms->resinfo[resind].chainid;
346
347 if (ch == 0)
348 {
349 ch = ' ';
350 }
351 }
352 if (resnr >= 10000)
353 {
354 resnr = resnr % 10000;
355 }
356 t_pdbinfo pdbinfo;
357 if (atoms->pdbinfo != nullptr)
358 {
359 pdbinfo = atoms->pdbinfo[i];
360 }
361 else
362 {
363 gmx_pdbinfo_init_default(&pdbinfo);
364 }
365 type = static_cast<enum PDB_record>(pdbinfo.type);
366 altloc = pdbinfo.altloc;
367 if (!isalnum(altloc))
368 {
369 altloc = ' ';
370 }
371 occup = bOccup ? 1.0 : pdbinfo.occup;
372 bfac = pdbinfo.bfac;
373
374 gmx_fprintf_pdb_atomline(out,
375 type,
376 i+1,
377 nm,
378 altloc,
379 resnm,
380 ch,
381 resnr,
382 resic,
383 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ],
384 occup,
385 bfac,
386 atoms->atom[i].elem);
387
388 if (atoms->pdbinfo && atoms->pdbinfo[i].bAnisotropic)
389 {
390 fprintf(out, "ANISOU%5d %-4.4s%4.4s%c%4d%c %7d%7d%7d%7d%7d%7d\n",
391 (i+1)%100000, nm, resnm, ch, resnr,
392 (resic == '\0') ? ' ' : resic,
393 atoms->pdbinfo[i].uij[0], atoms->pdbinfo[i].uij[1],
394 atoms->pdbinfo[i].uij[2], atoms->pdbinfo[i].uij[3],
395 atoms->pdbinfo[i].uij[4], atoms->pdbinfo[i].uij[5]);
396 }
397 }
398
399 fprintf(out, "TER\n");
400 fprintf(out, "ENDMDL\n");
401
402 if (NULL != gc)
403 {
404 /* Write conect records */
405 for (i = 0; (i < gc->nconect); i++)
406 {
407 fprintf(out, "CONECT%5d%5d\n", gc->conect[i].ai+1, gc->conect[i].aj+1);
408 }
409 }
410
411 gmx_residuetype_destroy(rt);
412 }
413
414 void write_pdbfile(FILE *out, const char *title, const t_atoms *atoms, const rvec x[],
415 int ePBC, const matrix box, char chainid, int model_nr, gmx_conect conect, gmx_bool bTerSepChains)
416 {
417 int i, *index;
418
419 snew(index, atoms->nr);
420 for (i = 0; i < atoms->nr; i++)
421 {
422 index[i] = i;
423 }
424 write_pdbfile_indexed(out, title, atoms, x, ePBC, box, chainid, model_nr,
425 atoms->nr, index, conect, bTerSepChains);
426 sfree(index);
427 }
428
429 int line2type(char *line)
430 {
431 int k;
432 char type[8];
433
434 for (k = 0; (k < 6); k++)
435 {
436 type[k] = line[k];
437 }
438 type[k] = '\0';
439
440 for (k = 0; (k < epdbNR); k++)
441 {
442 if (std::strncmp(type, pdbtp[k], strlen(pdbtp[k])) == 0)
443 {
444 break;
445 }
446 }
447
448 return k;
449 }
450
451 static void read_anisou(char line[], int natom, t_atoms *atoms)
452 {
453 int i, j, k, atomnr;
454 char nc = '\0';
455 char anr[12], anm[12];
456
457 /* Skip over type */
458 j = 6;
459 for (k = 0; (k < 5); k++, j++)
460 {
461 anr[k] = line[j];
462 }
463 anr[k] = nc;
464 j++;
465 for (k = 0; (k < 4); k++, j++)
466 {
467 anm[k] = line[j];
468 }
469 anm[k] = nc;
470 j++;
471
472 /* Strip off spaces */
473 trim(anm);
474
475 /* Search backwards for number and name only */
476 atomnr = std::strtol(anr, NULL, 10);
477 for (i = natom-1; (i >= 0); i--)
478 {
479 if ((std::strcmp(anm, *(atoms->atomname[i])) == 0) &&
480 (atomnr == atoms->pdbinfo[i].atomnr))
481 {
482 break;
483 }
484 }
485 if (i < 0)
486 {
487 fprintf(stderr, "Skipping ANISOU record (atom %s %d not found)\n",
488 anm, atomnr);
489 }
490 else
491 {
492 if (sscanf(line+29, "%d%d%d%d%d%d",
493 &atoms->pdbinfo[i].uij[U11], &atoms->pdbinfo[i].uij[U22],
494 &atoms->pdbinfo[i].uij[U33], &atoms->pdbinfo[i].uij[U12],
495 &atoms->pdbinfo[i].uij[U13], &atoms->pdbinfo[i].uij[U23])
496 == 6)
497 {
498 atoms->pdbinfo[i].bAnisotropic = TRUE;
499 }
500 else
501 {
502 fprintf(stderr, "Invalid ANISOU record for atom %d\n", i);
503 atoms->pdbinfo[i].bAnisotropic = FALSE;
504 }
505 }
506 }
507
508 void get_pdb_atomnumber(const t_atoms *atoms, gmx_atomprop_t aps)
509 {
510 int i, atomnumber, len;
511 size_t k;
512 char anm[6], anm_copy[6], *ptr;
513 char nc = '\0';
514 real eval;
515
516 if (!atoms->pdbinfo)
517 {
518 gmx_incons("Trying to deduce atomnumbers when no pdb information is present");
519 }
520 for (i = 0; (i < atoms->nr); i++)
521 {
522 std::strcpy(anm, atoms->pdbinfo[i].atomnm);
523 std::strcpy(anm_copy, atoms->pdbinfo[i].atomnm);
524 bool atomNumberSet = false;
525 len = strlen(anm);
526 if ((anm[0] != ' ') && ((len <= 2) || !std::isdigit(anm[2])))
527 {
528 anm_copy[2] = nc;
529 if (gmx_atomprop_query(aps, epropElement, "???", anm_copy, &eval))
530 {
531 atomnumber = std::round(eval);
532 atomNumberSet = true;
533 }
534 else
535 {
536 anm_copy[1] = nc;
537 if (gmx_atomprop_query(aps, epropElement, "???", anm_copy, &eval))
538 {
539 atomnumber = std::round(eval);
540 atomNumberSet = true;
541 }
542 }
543 }
544 if (!atomNumberSet)
545 {
546 k = 0;
547 while ((k < std::strlen(anm)) && (std::isspace(anm[k]) || std::isdigit(anm[k])))
548 {
549 k++;
550 }
551 anm_copy[0] = anm[k];
552 anm_copy[1] = nc;
553 if (gmx_atomprop_query(aps, epropElement, "???", anm_copy, &eval))
554 {
555 atomnumber = std::round(eval);
556 atomNumberSet = true;
557 }
558 }
559 if (atomNumberSet)
560 {
561 atoms->atom[i].atomnumber = atomnumber;
562 ptr = gmx_atomprop_element(aps, atomnumber);
563 if (debug)
564 {
565 fprintf(debug, "Atomnumber for atom '%s' is %d\n",
566 anm, atomnumber);
567 }
568 }
569 else
570 {
571 ptr = NULL;
572 }
573 std::strncpy(atoms->atom[i].elem, ptr == NULL ? "" : ptr, 4);
574 }
575 }
576
577 static int read_atom(t_symtab *symtab,
578 char line[], int type, int natom,
579 t_atoms *atoms, rvec x[], int chainnum, gmx_bool bChange)
580 {
581 t_atom *atomn;
582 int j, k;
583 char nc = '\0';
584 char anr[12], anm[12], anm_copy[12], altloc, resnm[12], rnr[12], elem[3];
585 char xc[12], yc[12], zc[12], occup[12], bfac[12];
586 unsigned char resic;
587 char chainid;
588 int resnr, atomnumber;
589
590 if (natom >= atoms->nr)
591 {
592 gmx_fatal(FARGS, "\nFound more atoms (%d) in pdb file than expected (%d)",
593 natom+1, atoms->nr);
594 }
595
596 /* Skip over type */
597 j = 6;
598 for (k = 0; (k < 5); k++, j++)
599 {
600 anr[k] = line[j];
601 }
602 anr[k] = nc;
603 trim(anr);
604 j++;
605 for (k = 0; (k < 4); k++, j++)
606 {
607 anm[k] = line[j];
608 }
609 anm[k] = nc;
610 std::strcpy(anm_copy, anm);
611 rtrim(anm_copy);
612 atomnumber = 0;
613 trim(anm);
614 altloc = line[j];
615 j++;
616 for (k = 0; (k < 4); k++, j++)
617 {
618 resnm[k] = line[j];
619 }
620 resnm[k] = nc;
621 trim(resnm);
622
623 chainid = line[j];
624 j++;
625
626 for (k = 0; (k < 4); k++, j++)
627 {
628 rnr[k] = line[j];
629 }
630 rnr[k] = nc;
631 trim(rnr);
632 resnr = std::strtol(rnr, NULL, 10);
633 resic = line[j];
634 j += 4;
635
636 /* X,Y,Z Coordinate */
637 for (k = 0; (k < 8); k++, j++)
638 {
639 xc[k] = line[j];
640 }
641 xc[k] = nc;
642 for (k = 0; (k < 8); k++, j++)
643 {
644 yc[k] = line[j];
645 }
646 yc[k] = nc;
647 for (k = 0; (k < 8); k++, j++)
648 {
649 zc[k] = line[j];
650 }
651 zc[k] = nc;
652
653 /* Weight */
654 for (k = 0; (k < 6); k++, j++)
655 {
656 occup[k] = line[j];
657 }
658 occup[k] = nc;
659
660 /* B-Factor */
661 for (k = 0; (k < 7); k++, j++)
662 {
663 bfac[k] = line[j];
664 }
665 bfac[k] = nc;
666
667 /* 10 blanks */
668 j += 10;
669
670 /* Element name */
671 for (k = 0; (k < 2); k++, j++)
672 {
673 elem[k] = line[j];
674 }
675 elem[k] = nc;
676 trim(elem);
677
678 if (atoms->atom)
679 {
680 atomn = &(atoms->atom[natom]);
681 if ((natom == 0) ||
682 atoms->resinfo[atoms->atom[natom-1].resind].nr != resnr ||
683 atoms->resinfo[atoms->atom[natom-1].resind].ic != resic ||
684 (strcmp(*atoms->resinfo[atoms->atom[natom-1].resind].name, resnm) != 0))
685 {
686 if (natom == 0)
687 {
688 atomn->resind = 0;
689 }
690 else
691 {
692 atomn->resind = atoms->atom[natom-1].resind + 1;
693 }
694 atoms->nres = atomn->resind + 1;
695 t_atoms_set_resinfo(atoms, natom, symtab, resnm, resnr, resic, chainnum, chainid);
696 }
697 else
698 {
699 atomn->resind = atoms->atom[natom-1].resind;
700 }
701 if (bChange)
702 {
703 xlate_atomname_pdb2gmx(anm);
704 }
705 atoms->atomname[natom] = put_symtab(symtab, anm);
706 atomn->m = 0.0;
707 atomn->q = 0.0;
708 atomn->atomnumber = atomnumber;
709 strncpy(atomn->elem, elem, 4);
710 }
711 x[natom][XX] = strtod(xc, NULL)*0.1;
712 x[natom][YY] = strtod(yc, NULL)*0.1;
713 x[natom][ZZ] = strtod(zc, NULL)*0.1;
714 if (atoms->pdbinfo)
715 {
716 atoms->pdbinfo[natom].type = type;
717 atoms->pdbinfo[natom].atomnr = strtol(anr, NULL, 10);
718 atoms->pdbinfo[natom].altloc = altloc;
719 strcpy(atoms->pdbinfo[natom].atomnm, anm_copy);
720 atoms->pdbinfo[natom].bfac = strtod(bfac, NULL);
721 atoms->pdbinfo[natom].occup = strtod(occup, NULL);
722 }
723 natom++;
724
725 return natom;
726 }
727
728 gmx_bool is_hydrogen(const char *nm)
729 {
730 char buf[30];
731
732 std::strcpy(buf, nm);
733 trim(buf);
734
735 if (buf[0] == 'H')
736 {
737 return TRUE;
738 }
739 else if ((std::isdigit(buf[0])) && (buf[1] == 'H'))
740 {
741 return TRUE;
742 }
743 return FALSE;
744 }
745
746 gmx_bool is_dummymass(const char *nm)
747 {
748 char buf[30];
749
750 std::strcpy(buf, nm);
751 trim(buf);
752
753 if ((buf[0] == 'M') && std::isdigit(buf[strlen(buf)-1]))
754 {
755 return TRUE;
756 }
757
758 return FALSE;
759 }
760
761 static void gmx_conect_addline(gmx_conect_t *con, char *line)
762 {
763 int n, ai, aj;
764 char format[32], form2[32];
765
766 sprintf(form2, "%%*s");
767 sprintf(format, "%s%%d", form2);
768 if (sscanf(line, format, &ai) == 1)
769 {
770 do
771 {
772 std::strcat(form2, "%*s");
773 sprintf(format, "%s%%d", form2);
774 n = sscanf(line, format, &aj);
775 if (n == 1)
776 {
777 srenew(con->conect, ++con->nconect);
778 con->conect[con->nconect-1].ai = ai-1;
779 con->conect[con->nconect-1].aj = aj-1;
780 }
781 }
782 while (n == 1);
783 }
784 }
785
786 void gmx_conect_dump(FILE *fp, gmx_conect conect)
787 {
788 gmx_conect_t *gc = (gmx_conect_t *)conect;
789 int i;
790
791 for (i = 0; (i < gc->nconect); i++)
792 {
793 fprintf(fp, "%6s%5d%5d\n", "CONECT",
794 gc->conect[i].ai+1, gc->conect[i].aj+1);
795 }
796 }
797
798 gmx_conect gmx_conect_init()
799 {
800 gmx_conect_t *gc;
801
802 snew(gc, 1);
803
804 return gc;
805 }
806
807 void gmx_conect_done(gmx_conect conect)
808 {
809 gmx_conect_t *gc = conect;
810
811 sfree(gc->conect);
812 }
813
814 gmx_bool gmx_conect_exist(gmx_conect conect, int ai, int aj)
815 {
816 gmx_conect_t *gc = conect;
817 int i;
818
819 /* if (!gc->bSorted)
820 sort_conect(gc);*/
821
822 for (i = 0; (i < gc->nconect); i++)
823 {
824 if (((gc->conect[i].ai == ai) &&
825 (gc->conect[i].aj == aj)) ||
826 ((gc->conect[i].aj == ai) &&
827 (gc->conect[i].ai == aj)))
828 {
829 return TRUE;
830 }
831 }
832 return FALSE;
833 }
834
835 void gmx_conect_add(gmx_conect conect, int ai, int aj)
836 {
837 gmx_conect_t *gc = (gmx_conect_t *)conect;
838
839 /* if (!gc->bSorted)
840 sort_conect(gc);*/
841
842 if (!gmx_conect_exist(conect, ai, aj))
843 {
844 srenew(gc->conect, ++gc->nconect);
845 gc->conect[gc->nconect-1].ai = ai;
846 gc->conect[gc->nconect-1].aj = aj;
847 }
848 }
849
850 int read_pdbfile(FILE *in, char *title, int *model_nr,
851 t_atoms *atoms, t_symtab *symtab, rvec x[], int *ePBC,
852 matrix box, gmx_bool bChange, gmx_conect conect)
853 {
854 gmx_conect_t *gc = conect;
855 gmx_bool bCOMPND;
856 gmx_bool bConnWarn = FALSE;
857 char line[STRLEN+1];
858 int line_type;
859 char *c, *d;
860 int natom, chainnum;
861 gmx_bool bStop = FALSE;
862
863 if (ePBC)
864 {
865 /* Only assume pbc when there is a CRYST1 entry */
866 *ePBC = epbcNONE;
867 }
868 if (box != NULL)
869 {
870 clear_mat(box);
871 }
872
873 atoms->haveMass = FALSE;
874 atoms->haveCharge = FALSE;
875 atoms->haveType = FALSE;
876 atoms->haveBState = FALSE;
877 atoms->havePdbInfo = (atoms->pdbinfo != NULL);
878
879 bCOMPND = FALSE;
880 title[0] = '\0';
881 natom = 0;
882 chainnum = 0;
883 while (!bStop && (fgets2(line, STRLEN, in) != NULL))
884 {
885 line_type = line2type(line);
886
887 switch (line_type)
888 {
889 case epdbATOM:
890 case epdbHETATM:
891 natom = read_atom(symtab, line, line_type, natom, atoms, x, chainnum, bChange);
892 break;
893
894 case epdbANISOU:
895 if (atoms->havePdbInfo)
896 {
897 read_anisou(line, natom, atoms);
898 }
899 break;
900
901 case epdbCRYST1:
902 read_cryst1(line, ePBC, box);
903 break;
904
905 case epdbTITLE:
906 case epdbHEADER:
907 if (std::strlen(line) > 6)
908 {
909 c = line+6;
910 /* skip HEADER or TITLE and spaces */
911 while (c[0] != ' ')
912 {
913 c++;
914 }
915 while (c[0] == ' ')
916 {
917 c++;
918 }
919 /* truncate after title */
920 d = std::strstr(c, " ");
921 if (d)
922 {
923 d[0] = '\0';
924 }
925 if (std::strlen(c) > 0)
926 {
927 std::strcpy(title, c);
928 }
929 }
930 break;
931
932 case epdbCOMPND:
933 if ((!std::strstr(line, ": ")) || (std::strstr(line+6, "MOLECULE:")))
934 {
935 if (!(c = std::strstr(line+6, "MOLECULE:")) )
936 {
937 c = line;
938 }
939 /* skip 'MOLECULE:' and spaces */
940 while (c[0] != ' ')
941 {
942 c++;
943 }
944 while (c[0] == ' ')
945 {
946 c++;
947 }
948 /* truncate after title */
949 d = strstr(c, " ");
950 if (d)
951 {
952 while ( (d[-1] == ';') && d > c)
953 {
954 d--;
955 }
956 d[0] = '\0';
957 }
958 if (strlen(c) > 0)
959 {
960 if (bCOMPND)
961 {
962 std::strcat(title, "; ");
963 std::strcat(title, c);
964 }
965 else
966 {
967 std::strcpy(title, c);
968 }
969 }
970 bCOMPND = TRUE;
971 }
972 break;
973
974 case epdbTER:
975 chainnum++;
976 break;
977
978 case epdbMODEL:
979 if (model_nr)
980 {
981 sscanf(line, "%*s%d", model_nr);
982 }
983 break;
984
985 case epdbENDMDL:
986 bStop = TRUE;
987 break;
988 case epdbCONECT:
989 if (gc)
990 {
991 gmx_conect_addline(gc, line);
992 }
993 else if (!bConnWarn)
994 {
995 fprintf(stderr, "WARNING: all CONECT records are ignored\n");
996 bConnWarn = TRUE;
997 }
998 break;
999
1000 default:
1001 break;
1002 }
1003 }
1004
1005 return natom;
1006 }
1007
1008 void get_pdb_coordnum(FILE *in, int *natoms)
1009 {
1010 char line[STRLEN];
1011
1012 *natoms = 0;
1013 while (fgets2(line, STRLEN, in))
1014 {
1015 if (std::strncmp(line, "ENDMDL", 6) == 0)
1016 {
1017 break;
1018 }
1019 if ((std::strncmp(line, "ATOM ", 6) == 0) || (std::strncmp(line, "HETATM", 6) == 0))
1020 {
1021 (*natoms)++;
1022 }
1023 }
1024 }
1025
1026 void gmx_pdb_read_conf(const char *infile,
1027 t_symtab *symtab, char ***name, t_atoms *atoms,
1028 rvec x[], int *ePBC, matrix box)
1029 {
1030 FILE *in = gmx_fio_fopen(infile, "r");
1031 char title[STRLEN];
1032 read_pdbfile(in, title, NULL, atoms, symtab, x, ePBC, box, TRUE, NULL);
1033 *name = put_symtab(symtab, title);
1034 gmx_fio_fclose(in);
1035 }
1036
1037 gmx_conect gmx_conect_generate(const t_topology *top)
1038 {
1039 int f, i;
1040 gmx_conect gc;
1041
1042 /* Fill the conect records */
1043 gc = gmx_conect_init();
1044
1045 for (f = 0; (f < F_NRE); f++)
1046 {
1047 if (IS_CHEMBOND(f))
1048 {
1049 for (i = 0; (i < top->idef.il[f].nr); i += interaction_function[f].nratoms+1)
1050 {
1051 gmx_conect_add(gc, top->idef.il[f].iatoms[i+1],
1052 top->idef.il[f].iatoms[i+2]);
1053 }
1054 }
1055 }
1056 return gc;
1057 }
1058
1059 int
1060 gmx_fprintf_pdb_atomline(FILE * fp,
1061 enum PDB_record record,
1062 int atom_seq_number,
1063 const char * atom_name,
1064 char alternate_location,
1065 const char * res_name,
1066 char chain_id,
1067 int res_seq_number,
1068 char res_insertion_code,
1069 real x,
1070 real y,
1071 real z,
1072 real occupancy,
1073 real b_factor,
1074 const char * element)
1075 {
1076 char tmp_atomname[6], tmp_resname[6];
1077 gmx_bool start_name_in_col13;
1078 int n;
1079
1080 if (record != epdbATOM && record != epdbHETATM)
1081 {
1082 gmx_fatal(FARGS, "Can only print PDB atom lines as ATOM or HETATM records");
1083 }
1084
1085 /* Format atom name */
1086 if (atom_name != NULL)
1087 {
1088 /* If the atom name is an element name with two chars, it should start already in column 13.
1089 * Otherwise it should start in column 14, unless the name length is 4 chars.
1090 */
1091 if ( (element != NULL) && (std::strlen(element) >= 2) && (gmx_strncasecmp(atom_name, element, 2) == 0) )
1092 {
1093 start_name_in_col13 = TRUE;
1094 }
1095 else
1096 {
1097 start_name_in_col13 = (std::strlen(atom_name) >= 4);
1098 }
1099 snprintf(tmp_atomname, sizeof(tmp_atomname), start_name_in_col13 ? "" : " ");
1100 std::strncat(tmp_atomname, atom_name, 4);
1101 tmp_atomname[5] = '\0';
1102 }
1103 else
1104 {
1105 tmp_atomname[0] = '\0';
1106 }
1107
1108 /* Format residue name */
1109 std::strncpy(tmp_resname, (res_name != NULL) ? res_name : "", 4);
1110 /* Make sure the string is terminated if strlen was > 4 */
1111 tmp_resname[4] = '\0';
1112 /* String is properly terminated, so now we can use strcat. By adding a
1113 * space we can write it right-justified, and if the original name was
1114 * three characters or less there will be a space added on the right side.
1115 */
1116 std::strcat(tmp_resname, " ");
1117
1118 /* Truncate integers so they fit */
1119 atom_seq_number = atom_seq_number % 100000;
1120 res_seq_number = res_seq_number % 10000;
1121
1122 n = fprintf(fp,
1123 "%-6s%5d %-4.4s%c%4.4s%c%4d%c %8.3f%8.3f%8.3f%6.2f%6.2f %2s\n",
1124 pdbtp[record],
1125 atom_seq_number,
1126 tmp_atomname,
1127 alternate_location,
1128 tmp_resname,
1129 chain_id,
1130 res_seq_number,
1131 res_insertion_code,
1132 x, y, z,
1133 occupancy,
1134 b_factor,
1135 (element != NULL) ? element : "");
1136
1137 return n;
1138 }
The ultimate molecular dynamics simulation package