src/gromacs/fileio/pdbio.cpp

   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,2015,2016, 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 "gmxpre.h"
  38
  39 #include "pdbio.h"
  40
  41 #include <cctype>
  42 #include <cmath>
  43 #include <cstdlib>
  44 #include <cstring>
  45
  46 #include <string>
  47
  48 #include "gromacs/fileio/gmxfio.h"
  49 #include "gromacs/math/units.h"
  50 #include "gromacs/math/vec.h"
  51 #include "gromacs/pbcutil/pbc.h"
  52 #include "gromacs/topology/atomprop.h"
  53 #include "gromacs/topology/ifunc.h"
  54 #include "gromacs/topology/residuetypes.h"
  55 #include "gromacs/topology/symtab.h"
  56 #include "gromacs/topology/topology.h"
  57 #include "gromacs/utility/coolstuff.h"
  58 #include "gromacs/utility/cstringutil.h"
  59 #include "gromacs/utility/fatalerror.h"
  60 #include "gromacs/utility/futil.h"
  61 #include "gromacs/utility/smalloc.h"
  62 #include "gromacs/utility/snprintf.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 #define REMARK_SIM_BOX "REMARK    THIS IS A SIMULATION BOX"
  82
  83 static void xlate_atomname_pdb2gmx(char *name)
  84 {
  85     int  i, length;
  86     char temp;
  87
  88     length = std::strlen(name);
  89     if (length > 3 && std::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 = std::strlen(name);
 106     if (length > 3 && std::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, const 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*gmx_angle(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*gmx_angle(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*gmx_angle(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 = std::cos(alpha*DEG2RAD);
 224             }
 225             else
 226             {
 227                 cosa = 0;
 228             }
 229             if (beta != 90.0)
 230             {
 231                 cosb = std::cos(beta*DEG2RAD);
 232             }
 233             else
 234             {
 235                 cosb = 0;
 236             }
 237             if (gamma != 90.0)
 238             {
 239                 cosg = std::cos(gamma*DEG2RAD);
 240                 sing = std::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] = std::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                            const t_atoms *atoms, const rvec x[],
 264                            int ePBC, const matrix box, char chainid,
 265                            int model_nr, int nindex, const int index[],
 266                            gmx_conect conect, gmx_bool bTerSepChains)
 267 {
 268     gmx_conect_t     *gc = (gmx_conect_t *)conect;
 269     char              resnm[6], nm[6];
 270     int               i, ii;
 271     int               resind, resnr;
 272     enum PDB_record   type;
 273     unsigned char     resic, ch;
 274     char              altloc;
 275     real              occup, bfac;
 276     gmx_bool          bOccup;
 277     int               chainnum, lastchainnum;
 278     gmx_residuetype_t*rt;
 279     const char       *p_restype;
 280     const char       *p_lastrestype;
 281
 282     gmx_residuetype_init(&rt);
 283
 284     fprintf(out, "TITLE     %s\n", (title && title[0]) ? title : gmx::bromacs().c_str());
 285     if (box && ( norm2(box[XX]) || norm2(box[YY]) || norm2(box[ZZ]) ) )
 286     {
 287         gmx_write_pdb_box(out, ePBC, box);
 288     }
 289     if (atoms->havePdbInfo)
 290     {
 291         /* Check whether any occupancies are set, in that case leave it as is,
 292          * otherwise set them all to one
 293          */
 294         bOccup = TRUE;
 295         for (ii = 0; (ii < nindex) && bOccup; ii++)
 296         {
 297             i      = index[ii];
 298             bOccup = bOccup && (atoms->pdbinfo[i].occup == 0.0);
 299         }
 300     }
 301     else
 302     {
 303         bOccup = FALSE;
 304     }
 305
 306     fprintf(out, "MODEL %8d\n", model_nr > 0 ? model_nr : 1);
 307
 308     lastchainnum      = -1;
 309     p_restype         = NULL;
 310
 311     for (ii = 0; ii < nindex; ii++)
 312     {
 313         i             = index[ii];
 314         resind        = atoms->atom[i].resind;
 315         chainnum      = atoms->resinfo[resind].chainnum;
 316         p_lastrestype = p_restype;
 317         gmx_residuetype_get_type(rt, *atoms->resinfo[resind].name, &p_restype);
 318
 319         /* Add a TER record if we changed chain, and if either the previous or this chain is protein/DNA/RNA. */
 320         if (bTerSepChains && ii > 0 && chainnum != lastchainnum)
 321         {
 322             /* Only add TER if the previous chain contained protein/DNA/RNA. */
 323             if (gmx_residuetype_is_protein(rt, p_lastrestype) || gmx_residuetype_is_dna(rt, p_lastrestype) || gmx_residuetype_is_rna(rt, p_lastrestype))
 324             {
 325                 fprintf(out, "TER\n");
 326             }
 327             lastchainnum    = chainnum;
 328         }
 329
 330         strncpy(resnm, *atoms->resinfo[resind].name, sizeof(resnm)-1);
 331         resnm[sizeof(resnm)-1] = 0;
 332         strncpy(nm, *atoms->atomname[i], sizeof(nm)-1);
 333         nm[sizeof(nm)-1] = 0;
 334
 335         /* rename HG12 to 2HG1, etc. */
 336         xlate_atomname_gmx2pdb(nm);
 337         resnr = atoms->resinfo[resind].nr;
 338         resic = atoms->resinfo[resind].ic;
 339         if (chainid != ' ')
 340         {
 341             ch = chainid;
 342         }
 343         else
 344         {
 345             ch = atoms->resinfo[resind].chainid;
 346
 347             if (ch == 0)
 348             {
 349                 ch = ' ';
 350             }
 351         }
 352         if (resnr >= 10000)
 353         {
 354             resnr = resnr % 10000;
 355         }
 356         t_pdbinfo pdbinfo;
 357         if (atoms->pdbinfo != nullptr)
 358         {
 359             pdbinfo = atoms->pdbinfo[i];
 360         }
 361         else
 362         {
 363             gmx_pdbinfo_init_default(&pdbinfo);
 364         }
 365         type   = static_cast<enum PDB_record>(pdbinfo.type);
 366         altloc = pdbinfo.altloc;
 367         if (!isalnum(altloc))
 368         {
 369             altloc = ' ';
 370         }
 371         occup = bOccup ? 1.0 : pdbinfo.occup;
 372         bfac  = pdbinfo.bfac;
 373
 374         gmx_fprintf_pdb_atomline(out,
 375                                  type,
 376                                  i+1,
 377                                  nm,
 378                                  altloc,
 379                                  resnm,
 380                                  ch,
 381                                  resnr,
 382                                  resic,
 383                                  10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ],
 384                                  occup,
 385                                  bfac,
 386                                  atoms->atom[i].elem);
 387
 388         if (atoms->pdbinfo && atoms->pdbinfo[i].bAnisotropic)
 389         {
 390             fprintf(out, "ANISOU%5d  %-4.4s%4.4s%c%4d%c %7d%7d%7d%7d%7d%7d\n",
 391                     (i+1)%100000, nm, resnm, ch, resnr,
 392                     (resic == '\0') ? ' ' : resic,
 393                     atoms->pdbinfo[i].uij[0], atoms->pdbinfo[i].uij[1],
 394                     atoms->pdbinfo[i].uij[2], atoms->pdbinfo[i].uij[3],
 395                     atoms->pdbinfo[i].uij[4], atoms->pdbinfo[i].uij[5]);
 396         }
 397     }
 398
 399     fprintf(out, "TER\n");
 400     fprintf(out, "ENDMDL\n");
 401
 402     if (NULL != gc)
 403     {
 404         /* Write conect records */
 405         for (i = 0; (i < gc->nconect); i++)
 406         {
 407             fprintf(out, "CONECT%5d%5d\n", gc->conect[i].ai+1, gc->conect[i].aj+1);
 408         }
 409     }
 410
 411     gmx_residuetype_destroy(rt);
 412 }
 413
 414 void write_pdbfile(FILE *out, const char *title, const t_atoms *atoms, const rvec x[],
 415                    int ePBC, const matrix box, char chainid, int model_nr, gmx_conect conect, gmx_bool bTerSepChains)
 416 {
 417     int i, *index;
 418
 419     snew(index, atoms->nr);
 420     for (i = 0; i < atoms->nr; i++)
 421     {
 422         index[i] = i;
 423     }
 424     write_pdbfile_indexed(out, title, atoms, x, ePBC, box, chainid, model_nr,
 425                           atoms->nr, index, conect, bTerSepChains);
 426     sfree(index);
 427 }
 428
 429 int line2type(char *line)
 430 {
 431     int  k;
 432     char type[8];
 433
 434     for (k = 0; (k < 6); k++)
 435     {
 436         type[k] = line[k];
 437     }
 438     type[k] = '\0';
 439
 440     for (k = 0; (k < epdbNR); k++)
 441     {
 442         if (std::strncmp(type, pdbtp[k], strlen(pdbtp[k])) == 0)
 443         {
 444             break;
 445         }
 446     }
 447
 448     return k;
 449 }
 450
 451 static void read_anisou(char line[], int natom, t_atoms *atoms)
 452 {
 453     int  i, j, k, atomnr;
 454     char nc = '\0';
 455     char anr[12], anm[12];
 456
 457     /* Skip over type */
 458     j = 6;
 459     for (k = 0; (k < 5); k++, j++)
 460     {
 461         anr[k] = line[j];
 462     }
 463     anr[k] = nc;
 464     j++;
 465     for (k = 0; (k < 4); k++, j++)
 466     {
 467         anm[k] = line[j];
 468     }
 469     anm[k] = nc;
 470     j++;
 471
 472     /* Strip off spaces */
 473     trim(anm);
 474
 475     /* Search backwards for number and name only */
 476     atomnr = std::strtol(anr, NULL, 10);
 477     for (i = natom-1; (i >= 0); i--)
 478     {
 479         if ((std::strcmp(anm, *(atoms->atomname[i])) == 0) &&
 480             (atomnr == atoms->pdbinfo[i].atomnr))
 481         {
 482             break;
 483         }
 484     }
 485     if (i < 0)
 486     {
 487         fprintf(stderr, "Skipping ANISOU record (atom %s %d not found)\n",
 488                 anm, atomnr);
 489     }
 490     else
 491     {
 492         if (sscanf(line+29, "%d%d%d%d%d%d",
 493                    &atoms->pdbinfo[i].uij[U11], &atoms->pdbinfo[i].uij[U22],
 494                    &atoms->pdbinfo[i].uij[U33], &atoms->pdbinfo[i].uij[U12],
 495                    &atoms->pdbinfo[i].uij[U13], &atoms->pdbinfo[i].uij[U23])
 496             == 6)
 497         {
 498             atoms->pdbinfo[i].bAnisotropic = TRUE;
 499         }
 500         else
 501         {
 502             fprintf(stderr, "Invalid ANISOU record for atom %d\n", i);
 503             atoms->pdbinfo[i].bAnisotropic = FALSE;
 504         }
 505     }
 506 }
 507
 508 void get_pdb_atomnumber(const t_atoms *atoms, gmx_atomprop_t aps)
 509 {
 510     int    i, atomnumber, len;
 511     size_t k;
 512     char   anm[6], anm_copy[6], *ptr;
 513     char   nc = '\0';
 514     real   eval;
 515
 516     if (!atoms->pdbinfo)
 517     {
 518         gmx_incons("Trying to deduce atomnumbers when no pdb information is present");
 519     }
 520     for (i = 0; (i < atoms->nr); i++)
 521     {
 522         std::strcpy(anm, atoms->pdbinfo[i].atomnm);
 523         std::strcpy(anm_copy, atoms->pdbinfo[i].atomnm);
 524         bool atomNumberSet = false;
 525         len        = strlen(anm);
 526         if ((anm[0] != ' ') && ((len <= 2) || !std::isdigit(anm[2])))
 527         {
 528             anm_copy[2] = nc;
 529             if (gmx_atomprop_query(aps, epropElement, "???", anm_copy, &eval))
 530             {
 531                 atomnumber    = std::round(eval);
 532                 atomNumberSet = true;
 533             }
 534             else
 535             {
 536                 anm_copy[1] = nc;
 537                 if (gmx_atomprop_query(aps, epropElement, "???", anm_copy, &eval))
 538                 {
 539                     atomnumber    = std::round(eval);
 540                     atomNumberSet = true;
 541                 }
 542             }
 543         }
 544         if (!atomNumberSet)
 545         {
 546             k = 0;
 547             while ((k < std::strlen(anm)) && (std::isspace(anm[k]) || std::isdigit(anm[k])))
 548             {
 549                 k++;
 550             }
 551             anm_copy[0] = anm[k];
 552             anm_copy[1] = nc;
 553             if (gmx_atomprop_query(aps, epropElement, "???", anm_copy, &eval))
 554             {
 555                 atomnumber    = std::round(eval);
 556                 atomNumberSet = true;
 557             }
 558         }
 559         if (atomNumberSet)
 560         {
 561             atoms->atom[i].atomnumber = atomnumber;
 562             ptr = gmx_atomprop_element(aps, atomnumber);
 563             if (debug)
 564             {
 565                 fprintf(debug, "Atomnumber for atom '%s' is %d\n",
 566                         anm, atomnumber);
 567             }
 568         }
 569         else
 570         {
 571             ptr = NULL;
 572         }
 573         std::strncpy(atoms->atom[i].elem, ptr == NULL ? "" : ptr, 4);
 574     }
 575 }
 576
 577 static int read_atom(t_symtab *symtab,
 578                      char line[], int type, int natom,
 579                      t_atoms *atoms, rvec x[], int chainnum, gmx_bool bChange)
 580 {
 581     t_atom       *atomn;
 582     int           j, k;
 583     char          nc = '\0';
 584     char          anr[12], anm[12], anm_copy[12], altloc, resnm[12], rnr[12], elem[3];
 585     char          xc[12], yc[12], zc[12], occup[12], bfac[12];
 586     unsigned char resic;
 587     char          chainid;
 588     int           resnr, atomnumber;
 589
 590     if (natom >= atoms->nr)
 591     {
 592         gmx_fatal(FARGS, "\nFound more atoms (%d) in pdb file than expected (%d)",
 593                   natom+1, atoms->nr);
 594     }
 595
 596     /* Skip over type */
 597     j = 6;
 598     for (k = 0; (k < 5); k++, j++)
 599     {
 600         anr[k] = line[j];
 601     }
 602     anr[k] = nc;
 603     trim(anr);
 604     j++;
 605     for (k = 0; (k < 4); k++, j++)
 606     {
 607         anm[k] = line[j];
 608     }
 609     anm[k] = nc;
 610     std::strcpy(anm_copy, anm);
 611     rtrim(anm_copy);
 612     atomnumber = 0;
 613     trim(anm);
 614     altloc = line[j];
 615     j++;
 616     for (k = 0; (k < 4); k++, j++)
 617     {
 618         resnm[k] = line[j];
 619     }
 620     resnm[k] = nc;
 621     trim(resnm);
 622
 623     chainid = line[j];
 624     j++;
 625
 626     for (k = 0; (k < 4); k++, j++)
 627     {
 628         rnr[k] = line[j];
 629     }
 630     rnr[k] = nc;
 631     trim(rnr);
 632     resnr = std::strtol(rnr, NULL, 10);
 633     resic = line[j];
 634     j    += 4;
 635
 636     /* X,Y,Z Coordinate */
 637     for (k = 0; (k < 8); k++, j++)
 638     {
 639         xc[k] = line[j];
 640     }
 641     xc[k] = nc;
 642     for (k = 0; (k < 8); k++, j++)
 643     {
 644         yc[k] = line[j];
 645     }
 646     yc[k] = nc;
 647     for (k = 0; (k < 8); k++, j++)
 648     {
 649         zc[k] = line[j];
 650     }
 651     zc[k] = nc;
 652
 653     /* Weight */
 654     for (k = 0; (k < 6); k++, j++)
 655     {
 656         occup[k] = line[j];
 657     }
 658     occup[k] = nc;
 659
 660     /* B-Factor */
 661     for (k = 0; (k < 7); k++, j++)
 662     {
 663         bfac[k] = line[j];
 664     }
 665     bfac[k] = nc;
 666
 667     /* 10 blanks */
 668     j += 10;
 669
 670     /* Element name */
 671     for (k = 0; (k < 2); k++, j++)
 672     {
 673         elem[k] = line[j];
 674     }
 675     elem[k] = nc;
 676     trim(elem);
 677
 678     if (atoms->atom)
 679     {
 680         atomn = &(atoms->atom[natom]);
 681         if ((natom == 0) ||
 682             atoms->resinfo[atoms->atom[natom-1].resind].nr != resnr ||
 683             atoms->resinfo[atoms->atom[natom-1].resind].ic != resic ||
 684             (strcmp(*atoms->resinfo[atoms->atom[natom-1].resind].name, resnm) != 0))
 685         {
 686             if (natom == 0)
 687             {
 688                 atomn->resind = 0;
 689             }
 690             else
 691             {
 692                 atomn->resind = atoms->atom[natom-1].resind + 1;
 693             }
 694             atoms->nres = atomn->resind + 1;
 695             t_atoms_set_resinfo(atoms, natom, symtab, resnm, resnr, resic, chainnum, chainid);
 696         }
 697         else
 698         {
 699             atomn->resind = atoms->atom[natom-1].resind;
 700         }
 701         if (bChange)
 702         {
 703             xlate_atomname_pdb2gmx(anm);
 704         }
 705         atoms->atomname[natom] = put_symtab(symtab, anm);
 706         atomn->m               = 0.0;
 707         atomn->q               = 0.0;
 708         atomn->atomnumber      = atomnumber;
 709         strncpy(atomn->elem, elem, 4);
 710     }
 711     x[natom][XX] = strtod(xc, NULL)*0.1;
 712     x[natom][YY] = strtod(yc, NULL)*0.1;
 713     x[natom][ZZ] = strtod(zc, NULL)*0.1;
 714     if (atoms->pdbinfo)
 715     {
 716         atoms->pdbinfo[natom].type   = type;
 717         atoms->pdbinfo[natom].atomnr = strtol(anr, NULL, 10);
 718         atoms->pdbinfo[natom].altloc = altloc;
 719         strcpy(atoms->pdbinfo[natom].atomnm, anm_copy);
 720         atoms->pdbinfo[natom].bfac  = strtod(bfac, NULL);
 721         atoms->pdbinfo[natom].occup = strtod(occup, NULL);
 722     }
 723     natom++;
 724
 725     return natom;
 726 }
 727
 728 gmx_bool is_hydrogen(const char *nm)
 729 {
 730     char buf[30];
 731
 732     std::strcpy(buf, nm);
 733     trim(buf);
 734
 735     if (buf[0] == 'H')
 736     {
 737         return TRUE;
 738     }
 739     else if ((std::isdigit(buf[0])) && (buf[1] == 'H'))
 740     {
 741         return TRUE;
 742     }
 743     return FALSE;
 744 }
 745
 746 gmx_bool is_dummymass(const char *nm)
 747 {
 748     char buf[30];
 749
 750     std::strcpy(buf, nm);
 751     trim(buf);
 752
 753     if ((buf[0] == 'M') && std::isdigit(buf[strlen(buf)-1]))
 754     {
 755         return TRUE;
 756     }
 757
 758     return FALSE;
 759 }
 760
 761 static void gmx_conect_addline(gmx_conect_t *con, char *line)
 762 {
 763     int  n, ai, aj;
 764     char format[32], form2[32];
 765
 766     sprintf(form2, "%%*s");
 767     sprintf(format, "%s%%d", form2);
 768     if (sscanf(line, format, &ai) == 1)
 769     {
 770         do
 771         {
 772             std::strcat(form2, "%*s");
 773             sprintf(format, "%s%%d", form2);
 774             n = sscanf(line, format, &aj);
 775             if (n == 1)
 776             {
 777                 srenew(con->conect, ++con->nconect);
 778                 con->conect[con->nconect-1].ai = ai-1;
 779                 con->conect[con->nconect-1].aj = aj-1;
 780             }
 781         }
 782         while (n == 1);
 783     }
 784 }
 785
 786 void gmx_conect_dump(FILE *fp, gmx_conect conect)
 787 {
 788     gmx_conect_t *gc = (gmx_conect_t *)conect;
 789     int           i;
 790
 791     for (i = 0; (i < gc->nconect); i++)
 792     {
 793         fprintf(fp, "%6s%5d%5d\n", "CONECT",
 794                 gc->conect[i].ai+1, gc->conect[i].aj+1);
 795     }
 796 }
 797
 798 gmx_conect gmx_conect_init()
 799 {
 800     gmx_conect_t *gc;
 801
 802     snew(gc, 1);
 803
 804     return gc;
 805 }
 806
 807 void gmx_conect_done(gmx_conect conect)
 808 {
 809     gmx_conect_t *gc = conect;
 810
 811     sfree(gc->conect);
 812 }
 813
 814 gmx_bool gmx_conect_exist(gmx_conect conect, int ai, int aj)
 815 {
 816     gmx_conect_t *gc = conect;
 817     int           i;
 818
 819     /* if (!gc->bSorted)
 820        sort_conect(gc);*/
 821
 822     for (i = 0; (i < gc->nconect); i++)
 823     {
 824         if (((gc->conect[i].ai == ai) &&
 825              (gc->conect[i].aj == aj)) ||
 826             ((gc->conect[i].aj == ai) &&
 827              (gc->conect[i].ai == aj)))
 828         {
 829             return TRUE;
 830         }
 831     }
 832     return FALSE;
 833 }
 834
 835 void gmx_conect_add(gmx_conect conect, int ai, int aj)
 836 {
 837     gmx_conect_t *gc = (gmx_conect_t *)conect;
 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, t_symtab *symtab, rvec x[], int *ePBC,
 852                  matrix box, gmx_bool bChange, gmx_conect conect)
 853 {
 854     gmx_conect_t *gc = conect;
 855     gmx_bool      bCOMPND;
 856     gmx_bool      bConnWarn = FALSE;
 857     char          line[STRLEN+1];
 858     int           line_type;
 859     char         *c, *d;
 860     int           natom, chainnum;
 861     gmx_bool      bStop   = FALSE;
 862
 863     if (ePBC)
 864     {
 865         /* Only assume pbc when there is a CRYST1 entry */
 866         *ePBC = epbcNONE;
 867     }
 868     if (box != NULL)
 869     {
 870         clear_mat(box);
 871     }
 872
 873     atoms->haveMass    = FALSE;
 874     atoms->haveCharge  = FALSE;
 875     atoms->haveType    = FALSE;
 876     atoms->haveBState  = FALSE;
 877     atoms->havePdbInfo = (atoms->pdbinfo != NULL);
 878
 879     bCOMPND  = FALSE;
 880     title[0] = '\0';
 881     natom    = 0;
 882     chainnum = 0;
 883     while (!bStop && (fgets2(line, STRLEN, in) != NULL))
 884     {
 885         line_type = line2type(line);
 886
 887         switch (line_type)
 888         {
 889             case epdbATOM:
 890             case epdbHETATM:
 891                 natom = read_atom(symtab, line, line_type, natom, atoms, x, chainnum, bChange);
 892                 break;
 893
 894             case epdbANISOU:
 895                 if (atoms->havePdbInfo)
 896                 {
 897                     read_anisou(line, natom, atoms);
 898                 }
 899                 break;
 900
 901             case epdbCRYST1:
 902                 read_cryst1(line, ePBC, box);
 903                 break;
 904
 905             case epdbTITLE:
 906             case epdbHEADER:
 907                 if (std::strlen(line) > 6)
 908                 {
 909                     c = line+6;
 910                     /* skip HEADER or TITLE and spaces */
 911                     while (c[0] != ' ')
 912                     {
 913                         c++;
 914                     }
 915                     while (c[0] == ' ')
 916                     {
 917                         c++;
 918                     }
 919                     /* truncate after title */
 920                     d = std::strstr(c, "      ");
 921                     if (d)
 922                     {
 923                         d[0] = '\0';
 924                     }
 925                     if (std::strlen(c) > 0)
 926                     {
 927                         std::strcpy(title, c);
 928                     }
 929                 }
 930                 break;
 931
 932             case epdbCOMPND:
 933                 if ((!std::strstr(line, ": ")) || (std::strstr(line+6, "MOLECULE:")))
 934                 {
 935                     if (!(c = std::strstr(line+6, "MOLECULE:")) )
 936                     {
 937                         c = line;
 938                     }
 939                     /* skip 'MOLECULE:' and spaces */
 940                     while (c[0] != ' ')
 941                     {
 942                         c++;
 943                     }
 944                     while (c[0] == ' ')
 945                     {
 946                         c++;
 947                     }
 948                     /* truncate after title */
 949                     d = strstr(c, "   ");
 950                     if (d)
 951                     {
 952                         while ( (d[-1] == ';') && d > c)
 953                         {
 954                             d--;
 955                         }
 956                         d[0] = '\0';
 957                     }
 958                     if (strlen(c) > 0)
 959                     {
 960                         if (bCOMPND)
 961                         {
 962                             std::strcat(title, "; ");
 963                             std::strcat(title, c);
 964                         }
 965                         else
 966                         {
 967                             std::strcpy(title, c);
 968                         }
 969                     }
 970                     bCOMPND = TRUE;
 971                 }
 972                 break;
 973
 974             case epdbTER:
 975                 chainnum++;
 976                 break;
 977
 978             case epdbMODEL:
 979                 if (model_nr)
 980                 {
 981                     sscanf(line, "%*s%d", model_nr);
 982                 }
 983                 break;
 984
 985             case epdbENDMDL:
 986                 bStop = TRUE;
 987                 break;
 988             case epdbCONECT:
 989                 if (gc)
 990                 {
 991                     gmx_conect_addline(gc, line);
 992                 }
 993                 else if (!bConnWarn)
 994                 {
 995                     fprintf(stderr, "WARNING: all CONECT records are ignored\n");
 996                     bConnWarn = TRUE;
 997                 }
 998                 break;
 999
1000             default:
1001                 break;
1002         }
1003     }
1004
1005     return natom;
1006 }
1007
1008 void get_pdb_coordnum(FILE *in, int *natoms)
1009 {
1010     char line[STRLEN];
1011
1012     *natoms = 0;
1013     while (fgets2(line, STRLEN, in))
1014     {
1015         if (std::strncmp(line, "ENDMDL", 6) == 0)
1016         {
1017             break;
1018         }
1019         if ((std::strncmp(line, "ATOM  ", 6) == 0) || (std::strncmp(line, "HETATM", 6) == 0))
1020         {
1021             (*natoms)++;
1022         }
1023     }
1024 }
1025
1026 void gmx_pdb_read_conf(const char *infile,
1027                        t_symtab *symtab, char ***name, t_atoms *atoms,
1028                        rvec x[], int *ePBC, matrix box)
1029 {
1030     FILE *in = gmx_fio_fopen(infile, "r");
1031     char  title[STRLEN];
1032     read_pdbfile(in, title, NULL, atoms, symtab, x, ePBC, box, TRUE, NULL);
1033     *name = put_symtab(symtab, title);
1034     gmx_fio_fclose(in);
1035 }
1036
1037 gmx_conect gmx_conect_generate(const t_topology *top)
1038 {
1039     int        f, i;
1040     gmx_conect gc;
1041
1042     /* Fill the conect records */
1043     gc  = gmx_conect_init();
1044
1045     for (f = 0; (f < F_NRE); f++)
1046     {
1047         if (IS_CHEMBOND(f))
1048         {
1049             for (i = 0; (i < top->idef.il[f].nr); i += interaction_function[f].nratoms+1)
1050             {
1051                 gmx_conect_add(gc, top->idef.il[f].iatoms[i+1],
1052                                top->idef.il[f].iatoms[i+2]);
1053             }
1054         }
1055     }
1056     return gc;
1057 }
1058
1059 int
1060 gmx_fprintf_pdb_atomline(FILE *            fp,
1061                          enum PDB_record   record,
1062                          int               atom_seq_number,
1063                          const char *      atom_name,
1064                          char              alternate_location,
1065                          const char *      res_name,
1066                          char              chain_id,
1067                          int               res_seq_number,
1068                          char              res_insertion_code,
1069                          real              x,
1070                          real              y,
1071                          real              z,
1072                          real              occupancy,
1073                          real              b_factor,
1074                          const char *      element)
1075 {
1076     char     tmp_atomname[6], tmp_resname[6];
1077     gmx_bool start_name_in_col13;
1078     int      n;
1079
1080     if (record != epdbATOM && record != epdbHETATM)
1081     {
1082         gmx_fatal(FARGS, "Can only print PDB atom lines as ATOM or HETATM records");
1083     }
1084
1085     /* Format atom name */
1086     if (atom_name != NULL)
1087     {
1088         /* If the atom name is an element name with two chars, it should start already in column 13.
1089          * Otherwise it should start in column 14, unless the name length is 4 chars.
1090          */
1091         if ( (element != NULL) && (std::strlen(element) >= 2) && (gmx_strncasecmp(atom_name, element, 2) == 0) )
1092         {
1093             start_name_in_col13 = TRUE;
1094         }
1095         else
1096         {
1097             start_name_in_col13 = (std::strlen(atom_name) >= 4);
1098         }
1099         snprintf(tmp_atomname, sizeof(tmp_atomname), start_name_in_col13 ? "" : " ");
1100         std::strncat(tmp_atomname, atom_name, 4);
1101         tmp_atomname[5] = '\0';
1102     }
1103     else
1104     {
1105         tmp_atomname[0] = '\0';
1106     }
1107
1108     /* Format residue name */
1109     std::strncpy(tmp_resname, (res_name != NULL) ? res_name : "", 4);
1110     /* Make sure the string is terminated if strlen was > 4 */
1111     tmp_resname[4] = '\0';
1112     /* String is properly terminated, so now we can use strcat. By adding a
1113      * space we can write it right-justified, and if the original name was
1114      * three characters or less there will be a space added on the right side.
1115      */
1116     std::strcat(tmp_resname, " ");
1117
1118     /* Truncate integers so they fit */
1119     atom_seq_number = atom_seq_number % 100000;
1120     res_seq_number  = res_seq_number % 10000;
1121
1122     n = fprintf(fp,
1123                 "%-6s%5d %-4.4s%c%4.4s%c%4d%c   %8.3f%8.3f%8.3f%6.2f%6.2f          %2s\n",
1124                 pdbtp[record],
1125                 atom_seq_number,
1126                 tmp_atomname,
1127                 alternate_location,
1128                 tmp_resname,
1129                 chain_id,
1130                 res_seq_number,
1131                 res_insertion_code,
1132                 x, y, z,
1133                 occupancy,
1134                 b_factor,
1135                 (element != NULL) ? element : "");
1136
1137     return n;
1138 }