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