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