src/gromacs/gmxlib/pdbio.c

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