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,2017,2018,2019, 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 #define REMARK_SIM_BOX "REMARK    THIS IS A SIMULATION BOX"
  81
  82 static void xlate_atomname_pdb2gmx(char *name)
  83 {
  84     int  i, length;
  85     char temp;
  86
  87     length = std::strlen(name);
  88     if (length > 3 && std::isdigit(name[0]))
  89     {
  90         temp = name[0];
  91         for (i = 1; i < length; i++)
  92         {
  93             name[i-1] = name[i];
  94         }
  95         name[length-1] = temp;
  96     }
  97 }
  98
  99 // Deliberately taking a copy of name to return it later
 100 static std::string xlate_atomname_gmx2pdb(std::string name)
 101 {
 102     size_t length = name.size();
 103     if (length > 3 && std::isdigit(name[length-1]))
 104     {
 105         char temp = name[length-1];
 106         for (size_t i = length-1; i > 0; --i)
 107         {
 108             name[i] = name[i-1];
 109         }
 110         name[0] = temp;
 111     }
 112     return name;
 113 }
 114
 115
 116 void gmx_write_pdb_box(FILE *out, int ePBC, const matrix box)
 117 {
 118     real alpha, beta, gamma;
 119
 120     if (ePBC == -1)
 121     {
 122         ePBC = guess_ePBC(box);
 123     }
 124
 125     if (ePBC == epbcNONE)
 126     {
 127         return;
 128     }
 129
 130     if (norm2(box[YY])*norm2(box[ZZ]) != 0)
 131     {
 132         alpha = RAD2DEG*gmx_angle(box[YY], box[ZZ]);
 133     }
 134     else
 135     {
 136         alpha = 90;
 137     }
 138     if (norm2(box[XX])*norm2(box[ZZ]) != 0)
 139     {
 140         beta  = RAD2DEG*gmx_angle(box[XX], box[ZZ]);
 141     }
 142     else
 143     {
 144         beta  = 90;
 145     }
 146     if (norm2(box[XX])*norm2(box[YY]) != 0)
 147     {
 148         gamma = RAD2DEG*gmx_angle(box[XX], box[YY]);
 149     }
 150     else
 151     {
 152         gamma = 90;
 153     }
 154     fprintf(out, "REMARK    THIS IS A SIMULATION BOX\n");
 155     if (ePBC != epbcSCREW)
 156     {
 157         fprintf(out, "CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f %-11s%4d\n",
 158                 10*norm(box[XX]), 10*norm(box[YY]), 10*norm(box[ZZ]),
 159                 alpha, beta, gamma, "P 1", 1);
 160     }
 161     else
 162     {
 163         /* Double the a-vector length and write the correct space group */
 164         fprintf(out, "CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f %-11s%4d\n",
 165                 20*norm(box[XX]), 10*norm(box[YY]), 10*norm(box[ZZ]),
 166                 alpha, beta, gamma, "P 21 1 1", 1);
 167
 168     }
 169 }
 170
 171 static void read_cryst1(char *line, int *ePBC, matrix box)
 172 {
 173 #define SG_SIZE 11
 174     char   sa[12], sb[12], sc[12], sg[SG_SIZE+1], ident;
 175     double fa, fb, fc, alpha, beta, gamma, cosa, cosb, cosg, sing;
 176     int    syma, symb, symc;
 177     int    ePBC_file;
 178
 179     sscanf(line, "%*s%s%s%s%lf%lf%lf", sa, sb, sc, &alpha, &beta, &gamma);
 180
 181     ePBC_file = -1;
 182     if (strlen(line) >= 55)
 183     {
 184         strncpy(sg, line+55, SG_SIZE);
 185         sg[SG_SIZE] = '\0';
 186         ident       = ' ';
 187         syma        = 0;
 188         symb        = 0;
 189         symc        = 0;
 190         sscanf(sg, "%c %d %d %d", &ident, &syma, &symb, &symc);
 191         if (ident == 'P' && syma ==  1 && symb <= 1 && symc <= 1)
 192         {
 193             fc        = strtod(sc, nullptr)*0.1;
 194             ePBC_file = (fc > 0 ? epbcXYZ : epbcXY);
 195         }
 196         if (ident == 'P' && syma == 21 && symb == 1 && symc == 1)
 197         {
 198             ePBC_file = epbcSCREW;
 199         }
 200     }
 201     if (ePBC)
 202     {
 203         *ePBC = ePBC_file;
 204     }
 205
 206     if (box)
 207     {
 208         fa = strtod(sa, nullptr)*0.1;
 209         fb = strtod(sb, nullptr)*0.1;
 210         fc = strtod(sc, nullptr)*0.1;
 211         if (ePBC_file == epbcSCREW)
 212         {
 213             fa *= 0.5;
 214         }
 215         clear_mat(box);
 216         box[XX][XX] = fa;
 217         if ((alpha != 90.0) || (beta != 90.0) || (gamma != 90.0))
 218         {
 219             if (alpha != 90.0)
 220             {
 221                 cosa = std::cos(alpha*DEG2RAD);
 222             }
 223             else
 224             {
 225                 cosa = 0;
 226             }
 227             if (beta != 90.0)
 228             {
 229                 cosb = std::cos(beta*DEG2RAD);
 230             }
 231             else
 232             {
 233                 cosb = 0;
 234             }
 235             if (gamma != 90.0)
 236             {
 237                 cosg = std::cos(gamma*DEG2RAD);
 238                 sing = std::sin(gamma*DEG2RAD);
 239             }
 240             else
 241             {
 242                 cosg = 0;
 243                 sing = 1;
 244             }
 245             box[YY][XX] = fb*cosg;
 246             box[YY][YY] = fb*sing;
 247             box[ZZ][XX] = fc*cosb;
 248             box[ZZ][YY] = fc*(cosa - cosb*cosg)/sing;
 249             box[ZZ][ZZ] = std::sqrt(fc*fc
 250                                     - box[ZZ][XX]*box[ZZ][XX] - box[ZZ][YY]*box[ZZ][YY]);
 251         }
 252         else
 253         {
 254             box[YY][YY] = fb;
 255             box[ZZ][ZZ] = fc;
 256         }
 257     }
 258 }
 259
 260 static int
 261 gmx_fprintf_pqr_atomline(FILE *            fp,
 262                          enum PDB_record   record,
 263                          int               atom_seq_number,
 264                          const char *      atom_name,
 265                          const char *      res_name,
 266                          char              chain_id,
 267                          int               res_seq_number,
 268                          real              x,
 269                          real              y,
 270                          real              z,
 271                          real              occupancy,
 272                          real              b_factor)
 273 {
 274     GMX_RELEASE_ASSERT(record == epdbATOM || record == epdbHETATM,
 275                        "Can only print PQR atom lines as ATOM or HETATM records");
 276
 277     /* Check atom name */
 278     GMX_RELEASE_ASSERT(atom_name != nullptr,
 279                        "Need atom information to print pqr");
 280
 281     /* Check residue name */
 282     GMX_RELEASE_ASSERT(res_name != nullptr,
 283                        "Need residue information to print pqr");
 284
 285     /* Truncate integers so they fit */
 286     atom_seq_number = atom_seq_number % 100000;
 287     res_seq_number  = res_seq_number % 10000;
 288
 289     int n = fprintf(fp,
 290                     "%-6s%5d %-4.4s%4.4s%c%4d %8.3f %8.3f %8.3f %6.2f %6.2f\n",
 291                     pdbtp[record],
 292                     atom_seq_number,
 293                     atom_name,
 294                     res_name,
 295                     chain_id,
 296                     res_seq_number,
 297                     x, y, z,
 298                     occupancy,
 299                     b_factor);
 300
 301     return n;
 302 }
 303
 304 void write_pdbfile_indexed(FILE *out, const char *title,
 305                            const t_atoms *atoms, const rvec x[],
 306                            int ePBC, const matrix box, char chainid,
 307                            int model_nr, int nindex, const int index[],
 308                            gmx_conect conect,
 309                            bool usePqrFormat)
 310 {
 311     gmx_conect_t   *gc = static_cast<gmx_conect_t *>(conect);
 312     enum PDB_record type;
 313     char            altloc;
 314     real            occup, bfac;
 315     gmx_bool        bOccup;
 316
 317
 318     fprintf(out, "TITLE     %s\n", (title && title[0]) ? title : gmx::bromacs().c_str());
 319     if (box && ( (norm2(box[XX]) != 0.0F) || (norm2(box[YY]) != 0.0F) || (norm2(box[ZZ]) != 0.0F) ) )
 320     {
 321         gmx_write_pdb_box(out, ePBC, box);
 322     }
 323     if (atoms->havePdbInfo)
 324     {
 325         /* Check whether any occupancies are set, in that case leave it as is,
 326          * otherwise set them all to one
 327          */
 328         bOccup = TRUE;
 329         for (int ii = 0; (ii < nindex) && bOccup; ii++)
 330         {
 331             int i      = index[ii];
 332             bOccup = bOccup && (atoms->pdbinfo[i].occup == 0.0);
 333         }
 334     }
 335     else
 336     {
 337         bOccup = FALSE;
 338     }
 339
 340     fprintf(out, "MODEL %8d\n", model_nr > 0 ? model_nr : 1);
 341
 342     ResidueType rt;
 343     for (int ii = 0; ii < nindex; ii++)
 344     {
 345         int         i             = index[ii];
 346         int         resind        = atoms->atom[i].resind;
 347         std::string resnm         = *atoms->resinfo[resind].name;
 348         std::string nm            = *atoms->atomname[i];
 349
 350         /* rename HG12 to 2HG1, etc. */
 351         nm    = xlate_atomname_gmx2pdb(nm);
 352         int           resnr = atoms->resinfo[resind].nr;
 353         unsigned char resic = atoms->resinfo[resind].ic;
 354         unsigned char ch;
 355         if (chainid != ' ')
 356         {
 357             ch = chainid;
 358         }
 359         else
 360         {
 361             ch = atoms->resinfo[resind].chainid;
 362
 363             if (ch == 0)
 364             {
 365                 ch = ' ';
 366             }
 367         }
 368         if (resnr >= 10000)
 369         {
 370             resnr = resnr % 10000;
 371         }
 372         t_pdbinfo pdbinfo;
 373         if (atoms->pdbinfo != nullptr)
 374         {
 375             pdbinfo = atoms->pdbinfo[i];
 376         }
 377         else
 378         {
 379             gmx_pdbinfo_init_default(&pdbinfo);
 380         }
 381         type   = static_cast<enum PDB_record>(pdbinfo.type);
 382         altloc = pdbinfo.altloc;
 383         if (!isalnum(altloc))
 384         {
 385             altloc = ' ';
 386         }
 387         occup = bOccup ? 1.0 : pdbinfo.occup;
 388         bfac  = pdbinfo.bfac;
 389         if (!usePqrFormat)
 390         {
 391             gmx_fprintf_pdb_atomline(out,
 392                                      type,
 393                                      i+1,
 394                                      nm.c_str(),
 395                                      altloc,
 396                                      resnm.c_str(),
 397                                      ch,
 398                                      resnr,
 399                                      resic,
 400                                      10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ],
 401                                      occup,
 402                                      bfac,
 403                                      atoms->atom[i].elem);
 404
 405             if (atoms->pdbinfo && atoms->pdbinfo[i].bAnisotropic)
 406             {
 407                 fprintf(out, "ANISOU%5d  %-4.4s%4.4s%c%4d%c %7d%7d%7d%7d%7d%7d\n",
 408                         (i+1)%100000, nm.c_str(), resnm.c_str(), ch, resnr,
 409                         (resic == '\0') ? ' ' : resic,
 410                         atoms->pdbinfo[i].uij[0], atoms->pdbinfo[i].uij[1],
 411                         atoms->pdbinfo[i].uij[2], atoms->pdbinfo[i].uij[3],
 412                         atoms->pdbinfo[i].uij[4], atoms->pdbinfo[i].uij[5]);
 413             }
 414         }
 415         else
 416         {
 417             gmx_fprintf_pqr_atomline(out,
 418                                      type,
 419                                      i+1,
 420                                      nm.c_str(),
 421                                      resnm.c_str(),
 422                                      ch,
 423                                      resnr,
 424                                      10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ],
 425                                      occup,
 426                                      bfac);
 427         }
 428     }
 429
 430     fprintf(out, "TER\n");
 431     fprintf(out, "ENDMDL\n");
 432
 433     if (nullptr != gc)
 434     {
 435         /* Write conect records */
 436         for (int i = 0; (i < gc->nconect); i++)
 437         {
 438             fprintf(out, "CONECT%5d%5d\n", gc->conect[i].ai+1, gc->conect[i].aj+1);
 439         }
 440     }
 441 }
 442
 443 void write_pdbfile(FILE *out, const char *title, const t_atoms *atoms, const rvec x[],
 444                    int ePBC, const matrix box, char chainid, int model_nr, gmx_conect conect)
 445 {
 446     int 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, false);
 455     sfree(index);
 456 }
 457
 458 static int line2type(const 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 (std::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 = std::strtol(anr, nullptr, 10);
 506     for (i = natom-1; (i >= 0); i--)
 507     {
 508         if ((std::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(const t_atoms *atoms, AtomProperties *aps)
 538 {
 539     int    i, atomnumber, len;
 540     size_t k;
 541     char   anm[6], anm_copy[6];
 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         std::strcpy(anm, atoms->pdbinfo[i].atomnm);
 552         std::strcpy(anm_copy, atoms->pdbinfo[i].atomnm);
 553         bool atomNumberSet = false;
 554         len        = strlen(anm);
 555         if ((anm[0] != ' ') && ((len <= 2) || !std::isdigit(anm[2])))
 556         {
 557             anm_copy[2] = nc;
 558             if (aps->setAtomProperty(epropElement, "???", anm_copy, &eval))
 559             {
 560                 atomnumber    = gmx::roundToInt(eval);
 561                 atomNumberSet = true;
 562             }
 563             else
 564             {
 565                 anm_copy[1] = nc;
 566                 if (aps->setAtomProperty(epropElement, "???", anm_copy, &eval))
 567                 {
 568                     atomnumber    = gmx::roundToInt(eval);
 569                     atomNumberSet = true;
 570                 }
 571             }
 572         }
 573         if (!atomNumberSet)
 574         {
 575             k = 0;
 576             while ((k < std::strlen(anm)) && (std::isspace(anm[k]) || std::isdigit(anm[k])))
 577             {
 578                 k++;
 579             }
 580             anm_copy[0] = anm[k];
 581             anm_copy[1] = nc;
 582             if (aps->setAtomProperty(epropElement, "???", anm_copy, &eval))
 583             {
 584                 atomnumber    = gmx::roundToInt(eval);
 585                 atomNumberSet = true;
 586             }
 587         }
 588         std::string buf;
 589         if (atomNumberSet)
 590         {
 591             atoms->atom[i].atomnumber = atomnumber;
 592             buf = aps->elementFromAtomNumber(atomnumber);
 593             if (debug)
 594             {
 595                 fprintf(debug, "Atomnumber for atom '%s' is %d\n",
 596                         anm, atomnumber);
 597             }
 598         }
 599         buf.resize(3);
 600         std::strncpy(atoms->atom[i].elem, buf.c_str(), 4);
 601     }
 602 }
 603
 604 static int read_atom(t_symtab *symtab,
 605                      const char line[], int type, int natom,
 606                      t_atoms *atoms, rvec x[], int chainnum, gmx_bool bChange)
 607 {
 608     t_atom       *atomn;
 609     int           j, k;
 610     char          nc = '\0';
 611     char          anr[12], anm[12], anm_copy[12], altloc, resnm[12], rnr[12], elem[3];
 612     char          xc[12], yc[12], zc[12], occup[12], bfac[12];
 613     unsigned char resic;
 614     char          chainid;
 615     int           resnr, atomnumber;
 616
 617     if (natom >= atoms->nr)
 618     {
 619         gmx_fatal(FARGS, "\nFound more atoms (%d) in pdb file than expected (%d)",
 620                   natom+1, atoms->nr);
 621     }
 622
 623     /* Skip over type */
 624     j = 6;
 625     for (k = 0; (k < 5); k++, j++)
 626     {
 627         anr[k] = line[j];
 628     }
 629     anr[k] = nc;
 630     trim(anr);
 631     j++;
 632     for (k = 0; (k < 4); k++, j++)
 633     {
 634         anm[k] = line[j];
 635     }
 636     anm[k] = nc;
 637     std::strcpy(anm_copy, anm);
 638     rtrim(anm_copy);
 639     atomnumber = 0;
 640     trim(anm);
 641     altloc = line[j];
 642     j++;
 643     for (k = 0; (k < 4); k++, j++)
 644     {
 645         resnm[k] = line[j];
 646     }
 647     resnm[k] = nc;
 648     trim(resnm);
 649
 650     chainid = line[j];
 651     j++;
 652
 653     for (k = 0; (k < 4); k++, j++)
 654     {
 655         rnr[k] = line[j];
 656     }
 657     rnr[k] = nc;
 658     trim(rnr);
 659     resnr = std::strtol(rnr, nullptr, 10);
 660     resic = line[j];
 661     j    += 4;
 662
 663     /* X,Y,Z Coordinate */
 664     for (k = 0; (k < 8); k++, j++)
 665     {
 666         xc[k] = line[j];
 667     }
 668     xc[k] = nc;
 669     for (k = 0; (k < 8); k++, j++)
 670     {
 671         yc[k] = line[j];
 672     }
 673     yc[k] = nc;
 674     for (k = 0; (k < 8); k++, j++)
 675     {
 676         zc[k] = line[j];
 677     }
 678     zc[k] = nc;
 679
 680     /* Weight */
 681     for (k = 0; (k < 6); k++, j++)
 682     {
 683         occup[k] = line[j];
 684     }
 685     occup[k] = nc;
 686
 687     /* B-Factor */
 688     for (k = 0; (k < 7); k++, j++)
 689     {
 690         bfac[k] = line[j];
 691     }
 692     bfac[k] = nc;
 693
 694     /* 10 blanks */
 695     j += 10;
 696
 697     /* Element name */
 698     for (k = 0; (k < 2); k++, j++)
 699     {
 700         elem[k] = line[j];
 701     }
 702     elem[k] = nc;
 703     trim(elem);
 704
 705     if (atoms->atom)
 706     {
 707         atomn = &(atoms->atom[natom]);
 708         if ((natom == 0) ||
 709             atoms->resinfo[atoms->atom[natom-1].resind].nr != resnr ||
 710             atoms->resinfo[atoms->atom[natom-1].resind].ic != resic ||
 711             (strcmp(*atoms->resinfo[atoms->atom[natom-1].resind].name, resnm) != 0))
 712         {
 713             if (natom == 0)
 714             {
 715                 atomn->resind = 0;
 716             }
 717             else
 718             {
 719                 atomn->resind = atoms->atom[natom-1].resind + 1;
 720             }
 721             atoms->nres = atomn->resind + 1;
 722             t_atoms_set_resinfo(atoms, natom, symtab, resnm, resnr, resic, chainnum, chainid);
 723         }
 724         else
 725         {
 726             atomn->resind = atoms->atom[natom-1].resind;
 727         }
 728         if (bChange)
 729         {
 730             xlate_atomname_pdb2gmx(anm);
 731         }
 732         atoms->atomname[natom] = put_symtab(symtab, anm);
 733         atomn->m               = 0.0;
 734         atomn->q               = 0.0;
 735         atomn->atomnumber      = atomnumber;
 736         strncpy(atomn->elem, elem, 4);
 737     }
 738     x[natom][XX] = strtod(xc, nullptr)*0.1;
 739     x[natom][YY] = strtod(yc, nullptr)*0.1;
 740     x[natom][ZZ] = strtod(zc, nullptr)*0.1;
 741     if (atoms->pdbinfo)
 742     {
 743         atoms->pdbinfo[natom].type   = type;
 744         atoms->pdbinfo[natom].atomnr = strtol(anr, nullptr, 10);
 745         atoms->pdbinfo[natom].altloc = altloc;
 746         strcpy(atoms->pdbinfo[natom].atomnm, anm_copy);
 747         atoms->pdbinfo[natom].bfac  = strtod(bfac, nullptr);
 748         atoms->pdbinfo[natom].occup = strtod(occup, nullptr);
 749     }
 750     natom++;
 751
 752     return natom;
 753 }
 754
 755 gmx_bool is_hydrogen(const char *nm)
 756 {
 757     char buf[30];
 758
 759     std::strcpy(buf, nm);
 760     trim(buf);
 761
 762     if (buf[0] == 'H')
 763     {
 764         return TRUE;
 765     }
 766     else if ((std::isdigit(buf[0])) && (buf[1] == 'H'))
 767     {
 768         return TRUE;
 769     }
 770     return FALSE;
 771 }
 772
 773 gmx_bool is_dummymass(const char *nm)
 774 {
 775     char buf[30];
 776
 777     std::strcpy(buf, nm);
 778     trim(buf);
 779
 780     return (buf[0] == 'M') && (std::isdigit(buf[strlen(buf)-1]) != 0);
 781 }
 782
 783 static void gmx_conect_addline(gmx_conect_t *con, char *line)
 784 {
 785     int         n, ai, aj;
 786
 787     std::string form2  = "%%*s";
 788     std::string format = form2 + "%%d";
 789     if (sscanf(line, format.c_str(), &ai) == 1)
 790     {
 791         do
 792         {
 793             form2 += "%*s";
 794             format = form2 + "%%d";
 795             n      = sscanf(line, format.c_str(), &aj);
 796             if (n == 1)
 797             {
 798                 srenew(con->conect, ++con->nconect);
 799                 con->conect[con->nconect-1].ai = ai-1;
 800                 con->conect[con->nconect-1].aj = aj-1;
 801             }
 802         }
 803         while (n == 1);
 804     }
 805 }
 806
 807 void gmx_conect_dump(FILE *fp, gmx_conect conect)
 808 {
 809     gmx_conect_t *gc = static_cast<gmx_conect_t *>(conect);
 810     int           i;
 811
 812     for (i = 0; (i < gc->nconect); i++)
 813     {
 814         fprintf(fp, "%6s%5d%5d\n", "CONECT",
 815                 gc->conect[i].ai+1, gc->conect[i].aj+1);
 816     }
 817 }
 818
 819 gmx_conect gmx_conect_init()
 820 {
 821     gmx_conect_t *gc;
 822
 823     snew(gc, 1);
 824
 825     return gc;
 826 }
 827
 828 void gmx_conect_done(gmx_conect conect)
 829 {
 830     gmx_conect_t *gc = conect;
 831
 832     sfree(gc->conect);
 833 }
 834
 835 gmx_bool gmx_conect_exist(gmx_conect conect, int ai, int aj)
 836 {
 837     gmx_conect_t *gc = conect;
 838     int           i;
 839
 840     /* if (!gc->bSorted)
 841        sort_conect(gc);*/
 842
 843     for (i = 0; (i < gc->nconect); i++)
 844     {
 845         if (((gc->conect[i].ai == ai) &&
 846              (gc->conect[i].aj == aj)) ||
 847             ((gc->conect[i].aj == ai) &&
 848              (gc->conect[i].ai == aj)))
 849         {
 850             return TRUE;
 851         }
 852     }
 853     return FALSE;
 854 }
 855
 856 void gmx_conect_add(gmx_conect conect, int ai, int aj)
 857 {
 858     gmx_conect_t *gc = static_cast<gmx_conect_t *>(conect);
 859
 860     /* if (!gc->bSorted)
 861        sort_conect(gc);*/
 862
 863     if (!gmx_conect_exist(conect, ai, aj))
 864     {
 865         srenew(gc->conect, ++gc->nconect);
 866         gc->conect[gc->nconect-1].ai = ai;
 867         gc->conect[gc->nconect-1].aj = aj;
 868     }
 869 }
 870
 871 int read_pdbfile(FILE *in, char *title, int *model_nr,
 872                  t_atoms *atoms, t_symtab *symtab, rvec x[], int *ePBC,
 873                  matrix box, gmx_bool bChange, gmx_conect conect)
 874 {
 875     gmx_conect_t *gc = conect;
 876     gmx_bool      bCOMPND;
 877     gmx_bool      bConnWarn = FALSE;
 878     char          line[STRLEN+1];
 879     int           line_type;
 880     char         *c, *d;
 881     int           natom, chainnum;
 882     gmx_bool      bStop   = FALSE;
 883
 884     if (ePBC)
 885     {
 886         /* Only assume pbc when there is a CRYST1 entry */
 887         *ePBC = epbcNONE;
 888     }
 889     if (box != nullptr)
 890     {
 891         clear_mat(box);
 892     }
 893
 894     atoms->haveMass    = FALSE;
 895     atoms->haveCharge  = FALSE;
 896     atoms->haveType    = FALSE;
 897     atoms->haveBState  = FALSE;
 898     atoms->havePdbInfo = (atoms->pdbinfo != nullptr);
 899
 900     bCOMPND  = FALSE;
 901     title[0] = '\0';
 902     natom    = 0;
 903     chainnum = 0;
 904     while (!bStop && (fgets2(line, STRLEN, in) != nullptr))
 905     {
 906         line_type = line2type(line);
 907
 908         switch (line_type)
 909         {
 910             case epdbATOM:
 911             case epdbHETATM:
 912                 natom = read_atom(symtab, line, line_type, natom, atoms, x, chainnum, bChange);
 913                 break;
 914
 915             case epdbANISOU:
 916                 if (atoms->havePdbInfo)
 917                 {
 918                     read_anisou(line, natom, atoms);
 919                 }
 920                 break;
 921
 922             case epdbCRYST1:
 923                 read_cryst1(line, ePBC, box);
 924                 break;
 925
 926             case epdbTITLE:
 927             case epdbHEADER:
 928                 if (std::strlen(line) > 6)
 929                 {
 930                     c = line+6;
 931                     /* skip HEADER or TITLE and spaces */
 932                     while (c[0] != ' ')
 933                     {
 934                         c++;
 935                     }
 936                     while (c[0] == ' ')
 937                     {
 938                         c++;
 939                     }
 940                     /* truncate after title */
 941                     d = std::strstr(c, "      ");
 942                     if (d)
 943                     {
 944                         d[0] = '\0';
 945                     }
 946                     if (std::strlen(c) > 0)
 947                     {
 948                         std::strcpy(title, c);
 949                     }
 950                 }
 951                 break;
 952
 953             case epdbCOMPND:
 954                 if ((!std::strstr(line, ": ")) || (std::strstr(line+6, "MOLECULE:")))
 955                 {
 956                     if (!(c = std::strstr(line+6, "MOLECULE:")) )
 957                     {
 958                         c = line;
 959                     }
 960                     /* skip 'MOLECULE:' and spaces */
 961                     while (c[0] != ' ')
 962                     {
 963                         c++;
 964                     }
 965                     while (c[0] == ' ')
 966                     {
 967                         c++;
 968                     }
 969                     /* truncate after title */
 970                     d = strstr(c, "   ");
 971                     if (d)
 972                     {
 973                         while ( (d[-1] == ';') && d > c)
 974                         {
 975                             d--;
 976                         }
 977                         d[0] = '\0';
 978                     }
 979                     if (strlen(c) > 0)
 980                     {
 981                         if (bCOMPND)
 982                         {
 983                             std::strcat(title, "; ");
 984                             std::strcat(title, c);
 985                         }
 986                         else
 987                         {
 988                             std::strcpy(title, c);
 989                         }
 990                     }
 991                     bCOMPND = TRUE;
 992                 }
 993                 break;
 994
 995             case epdbTER:
 996                 chainnum++;
 997                 break;
 998
 999             case epdbMODEL:
1000                 if (model_nr)
1001                 {
1002                     sscanf(line, "%*s%d", model_nr);
1003                 }
1004                 break;
1005
1006             case epdbENDMDL:
1007                 bStop = TRUE;
1008                 break;
1009             case epdbCONECT:
1010                 if (gc)
1011                 {
1012                     gmx_conect_addline(gc, line);
1013                 }
1014                 else if (!bConnWarn)
1015                 {
1016                     fprintf(stderr, "WARNING: all CONECT records are ignored\n");
1017                     bConnWarn = TRUE;
1018                 }
1019                 break;
1020
1021             default:
1022                 break;
1023         }
1024     }
1025
1026     return natom;
1027 }
1028
1029 void get_pdb_coordnum(FILE *in, int *natoms)
1030 {
1031     char line[STRLEN];
1032
1033     *natoms = 0;
1034     while (fgets2(line, STRLEN, in))
1035     {
1036         if (std::strncmp(line, "ENDMDL", 6) == 0)
1037         {
1038             break;
1039         }
1040         if ((std::strncmp(line, "ATOM  ", 6) == 0) || (std::strncmp(line, "HETATM", 6) == 0))
1041         {
1042             (*natoms)++;
1043         }
1044     }
1045 }
1046
1047 void gmx_pdb_read_conf(const char *infile,
1048                        t_symtab *symtab, char **name, t_atoms *atoms,
1049                        rvec x[], int *ePBC, matrix box)
1050 {
1051     FILE *in = gmx_fio_fopen(infile, "r");
1052     char  title[STRLEN];
1053     read_pdbfile(in, title, nullptr, atoms, symtab, x, ePBC, box, TRUE, nullptr);
1054     if (name != nullptr)
1055     {
1056         *name = gmx_strdup(title);
1057     }
1058     gmx_fio_fclose(in);
1059 }
1060
1061 gmx_conect gmx_conect_generate(const t_topology *top)
1062 {
1063     int        f, i;
1064     gmx_conect gc;
1065
1066     /* Fill the conect records */
1067     gc  = gmx_conect_init();
1068
1069     for (f = 0; (f < F_NRE); f++)
1070     {
1071         if (IS_CHEMBOND(f))
1072         {
1073             for (i = 0; (i < top->idef.il[f].nr); i += interaction_function[f].nratoms+1)
1074             {
1075                 gmx_conect_add(gc, top->idef.il[f].iatoms[i+1],
1076                                top->idef.il[f].iatoms[i+2]);
1077             }
1078         }
1079     }
1080     return gc;
1081 }
1082
1083 int
1084 gmx_fprintf_pdb_atomline(FILE *            fp,
1085                          enum PDB_record   record,
1086                          int               atom_seq_number,
1087                          const char *      atom_name,
1088                          char              alternate_location,
1089                          const char *      res_name,
1090                          char              chain_id,
1091                          int               res_seq_number,
1092                          char              res_insertion_code,
1093                          real              x,
1094                          real              y,
1095                          real              z,
1096                          real              occupancy,
1097                          real              b_factor,
1098                          const char *      element)
1099 {
1100     char     tmp_atomname[6], tmp_resname[6];
1101     gmx_bool start_name_in_col13;
1102     int      n;
1103
1104     if (record != epdbATOM && record != epdbHETATM)
1105     {
1106         gmx_fatal(FARGS, "Can only print PDB atom lines as ATOM or HETATM records");
1107     }
1108
1109     /* Format atom name */
1110     if (atom_name != nullptr)
1111     {
1112         /* If the atom name is an element name with two chars, it should start already in column 13.
1113          * Otherwise it should start in column 14, unless the name length is 4 chars.
1114          */
1115         if ( (element != nullptr) && (std::strlen(element) >= 2) && (gmx_strncasecmp(atom_name, element, 2) == 0) )
1116         {
1117             start_name_in_col13 = TRUE;
1118         }
1119         else
1120         {
1121             start_name_in_col13 = (std::strlen(atom_name) >= 4);
1122         }
1123         snprintf(tmp_atomname, sizeof(tmp_atomname), start_name_in_col13 ? "" : " ");
1124         std::strncat(tmp_atomname, atom_name, 4);
1125         tmp_atomname[5] = '\0';
1126     }
1127     else
1128     {
1129         tmp_atomname[0] = '\0';
1130     }
1131
1132     /* Format residue name */
1133     std::strncpy(tmp_resname, (res_name != nullptr) ? res_name : "", 4);
1134     /* Make sure the string is terminated if strlen was > 4 */
1135     tmp_resname[4] = '\0';
1136     /* String is properly terminated, so now we can use strcat. By adding a
1137      * space we can write it right-justified, and if the original name was
1138      * three characters or less there will be a space added on the right side.
1139      */
1140     std::strcat(tmp_resname, " ");
1141
1142     /* Truncate integers so they fit */
1143     atom_seq_number = atom_seq_number % 100000;
1144     res_seq_number  = res_seq_number % 10000;
1145
1146     n = fprintf(fp,
1147                 "%-6s%5d %-4.4s%c%4.4s%c%4d%c   %8.3f%8.3f%8.3f%6.2f%6.2f          %2s\n",
1148                 pdbtp[record],
1149                 atom_seq_number,
1150                 tmp_atomname,
1151                 alternate_location,
1152                 tmp_resname,
1153                 chain_id,
1154                 res_seq_number,
1155                 res_insertion_code,
1156                 x, y, z,
1157                 occupancy,
1158                 b_factor,
1159                 (element != nullptr) ? element : "");
1160
1161     return n;
1162 }