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