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