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