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