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 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include "pdbio.h"
  41
  42 #include <cctype>
  43 #include <cmath>
  44 #include <cstdlib>
  45 #include <cstring>
  46
  47 #include <string>
  48
  49 #include "gromacs/fileio/gmxfio.h"
  50 #include "gromacs/math/units.h"
  51 #include "gromacs/math/vec.h"
  52 #include "gromacs/pbcutil/pbc.h"
  53 #include "gromacs/topology/atomprop.h"
  54 #include "gromacs/topology/ifunc.h"
  55 #include "gromacs/topology/residuetypes.h"
  56 #include "gromacs/topology/symtab.h"
  57 #include "gromacs/topology/topology.h"
  58 #include "gromacs/utility/coolstuff.h"
  59 #include "gromacs/utility/cstringutil.h"
  60 #include "gromacs/utility/fatalerror.h"
  61 #include "gromacs/utility/futil.h"
  62 #include "gromacs/utility/smalloc.h"
  63 #include "gromacs/utility/snprintf.h"
  64
  65 typedef struct
  66 {
  67     int ai, aj;
  68 } gmx_conection_t;
  69
  70 typedef struct gmx_conect_t
  71 {
  72     int              nconect;
  73     gmx_bool         bSorted;
  74     gmx_conection_t* conect;
  75 } gmx_conect_t;
  76
  77 static const char* pdbtp[epdbNR] = { "ATOM  ", "HETATM", "ANISOU", "CRYST1", "COMPND", "MODEL",
  78                                      "ENDMDL", "TER",    "HEADER", "TITLE",  "REMARK", "CONECT" };
  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", 10 * norm(box[XX]),
 158                 10 * norm(box[YY]), 10 * norm(box[ZZ]), alpha, beta, gamma, "P 1", 1);
 159     }
 160     else
 161     {
 162         /* Double the a-vector length and write the correct space group */
 163         fprintf(out, "CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f %-11s%4d\n", 20 * norm(box[XX]),
 164                 10 * norm(box[YY]), 10 * norm(box[ZZ]), alpha, beta, gamma, "P 21 1 1", 1);
 165     }
 166 }
 167
 168 static void read_cryst1(char* line, int* ePBC, matrix box)
 169 {
 170 #define SG_SIZE 11
 171     char   sa[12], sb[12], sc[12], sg[SG_SIZE + 1], ident;
 172     double fa, fb, fc, alpha, beta, gamma, cosa, cosb, cosg, sing;
 173     int    syma, symb, symc;
 174     int    ePBC_file;
 175
 176     sscanf(line, "%*s%s%s%s%lf%lf%lf", sa, sb, sc, &alpha, &beta, &gamma);
 177
 178     ePBC_file = -1;
 179     if (strlen(line) >= 55)
 180     {
 181         strncpy(sg, line + 55, SG_SIZE);
 182         sg[SG_SIZE] = '\0';
 183         ident       = ' ';
 184         syma        = 0;
 185         symb        = 0;
 186         symc        = 0;
 187         sscanf(sg, "%c %d %d %d", &ident, &syma, &symb, &symc);
 188         if (ident == 'P' && syma == 1 && symb <= 1 && symc <= 1)
 189         {
 190             fc        = strtod(sc, nullptr) * 0.1;
 191             ePBC_file = (fc > 0 ? epbcXYZ : epbcXY);
 192         }
 193         if (ident == 'P' && syma == 21 && symb == 1 && symc == 1)
 194         {
 195             ePBC_file = epbcSCREW;
 196         }
 197     }
 198     if (ePBC)
 199     {
 200         *ePBC = ePBC_file;
 201     }
 202
 203     if (box)
 204     {
 205         fa = strtod(sa, nullptr) * 0.1;
 206         fb = strtod(sb, nullptr) * 0.1;
 207         fc = strtod(sc, nullptr) * 0.1;
 208         if (ePBC_file == epbcSCREW)
 209         {
 210             fa *= 0.5;
 211         }
 212         clear_mat(box);
 213         box[XX][XX] = fa;
 214         if ((alpha != 90.0) || (beta != 90.0) || (gamma != 90.0))
 215         {
 216             if (alpha != 90.0)
 217             {
 218                 cosa = std::cos(alpha * DEG2RAD);
 219             }
 220             else
 221             {
 222                 cosa = 0;
 223             }
 224             if (beta != 90.0)
 225             {
 226                 cosb = std::cos(beta * DEG2RAD);
 227             }
 228             else
 229             {
 230                 cosb = 0;
 231             }
 232             if (gamma != 90.0)
 233             {
 234                 cosg = std::cos(gamma * DEG2RAD);
 235                 sing = std::sin(gamma * DEG2RAD);
 236             }
 237             else
 238             {
 239                 cosg = 0;
 240                 sing = 1;
 241             }
 242             box[YY][XX] = fb * cosg;
 243             box[YY][YY] = fb * sing;
 244             box[ZZ][XX] = fc * cosb;
 245             box[ZZ][YY] = fc * (cosa - cosb * cosg) / sing;
 246             box[ZZ][ZZ] = std::sqrt(fc * fc - box[ZZ][XX] * box[ZZ][XX] - box[ZZ][YY] * box[ZZ][YY]);
 247         }
 248         else
 249         {
 250             box[YY][YY] = fb;
 251             box[ZZ][ZZ] = fc;
 252         }
 253     }
 254 }
 255
 256 static int gmx_fprintf_pqr_atomline(FILE*           fp,
 257                                     enum PDB_record record,
 258                                     int             atom_seq_number,
 259                                     const char*     atom_name,
 260                                     const char*     res_name,
 261                                     char            chain_id,
 262                                     int             res_seq_number,
 263                                     real            x,
 264                                     real            y,
 265                                     real            z,
 266                                     real            occupancy,
 267                                     real            b_factor)
 268 {
 269     GMX_RELEASE_ASSERT(record == epdbATOM || record == epdbHETATM,
 270                        "Can only print PQR atom lines as ATOM or HETATM records");
 271
 272     /* Check atom name */
 273     GMX_RELEASE_ASSERT(atom_name != nullptr, "Need atom information to print pqr");
 274
 275     /* Check residue name */
 276     GMX_RELEASE_ASSERT(res_name != nullptr, "Need residue information to print pqr");
 277
 278     /* Truncate integers so they fit */
 279     atom_seq_number = atom_seq_number % 100000;
 280     res_seq_number  = res_seq_number % 10000;
 281
 282     int n = fprintf(fp, "%-6s%5d %-4.4s%4.4s%c%4d %8.3f %8.3f %8.3f %6.2f %6.2f\n", pdbtp[record],
 283                     atom_seq_number, atom_name, res_name, chain_id, res_seq_number, x, y, z,
 284                     occupancy, b_factor);
 285
 286     return n;
 287 }
 288
 289 void write_pdbfile_indexed(FILE*          out,
 290                            const char*    title,
 291                            const t_atoms* atoms,
 292                            const rvec     x[],
 293                            int            ePBC,
 294                            const matrix   box,
 295                            char           chainid,
 296                            int            model_nr,
 297                            int            nindex,
 298                            const int      index[],
 299                            gmx_conect     conect,
 300                            bool           usePqrFormat)
 301 {
 302     gmx_conect_t*   gc = static_cast<gmx_conect_t*>(conect);
 303     enum PDB_record type;
 304     char            altloc;
 305     real            occup, bfac;
 306     gmx_bool        bOccup;
 307
 308
 309     fprintf(out, "TITLE     %s\n", (title && title[0]) ? title : gmx::bromacs().c_str());
 310     if (box && ((norm2(box[XX]) != 0.0F) || (norm2(box[YY]) != 0.0F) || (norm2(box[ZZ]) != 0.0F)))
 311     {
 312         gmx_write_pdb_box(out, ePBC, box);
 313     }
 314     if (atoms->havePdbInfo)
 315     {
 316         /* Check whether any occupancies are set, in that case leave it as is,
 317          * otherwise set them all to one
 318          */
 319         bOccup = TRUE;
 320         for (int ii = 0; (ii < nindex) && bOccup; ii++)
 321         {
 322             int i  = index[ii];
 323             bOccup = bOccup && (atoms->pdbinfo[i].occup == 0.0);
 324         }
 325     }
 326     else
 327     {
 328         bOccup = FALSE;
 329     }
 330
 331     fprintf(out, "MODEL %8d\n", model_nr > 0 ? model_nr : 1);
 332
 333     ResidueType rt;
 334     for (int ii = 0; ii < nindex; ii++)
 335     {
 336         int         i      = index[ii];
 337         int         resind = atoms->atom[i].resind;
 338         std::string resnm  = *atoms->resinfo[resind].name;
 339         std::string nm     = *atoms->atomname[i];
 340
 341         /* rename HG12 to 2HG1, etc. */
 342         nm                  = xlate_atomname_gmx2pdb(nm);
 343         int           resnr = atoms->resinfo[resind].nr;
 344         unsigned char resic = atoms->resinfo[resind].ic;
 345         unsigned char ch;
 346         if (chainid != ' ')
 347         {
 348             ch = chainid;
 349         }
 350         else
 351         {
 352             ch = atoms->resinfo[resind].chainid;
 353
 354             if (ch == 0)
 355             {
 356                 ch = ' ';
 357             }
 358         }
 359         if (resnr >= 10000)
 360         {
 361             resnr = resnr % 10000;
 362         }
 363         t_pdbinfo pdbinfo;
 364         if (atoms->pdbinfo != nullptr)
 365         {
 366             pdbinfo = atoms->pdbinfo[i];
 367         }
 368         else
 369         {
 370             gmx_pdbinfo_init_default(&pdbinfo);
 371         }
 372         type   = static_cast<enum PDB_record>(pdbinfo.type);
 373         altloc = pdbinfo.altloc;
 374         if (!isalnum(altloc))
 375         {
 376             altloc = ' ';
 377         }
 378         occup = bOccup ? 1.0 : pdbinfo.occup;
 379         bfac  = pdbinfo.bfac;
 380         if (!usePqrFormat)
 381         {
 382             gmx_fprintf_pdb_atomline(out, type, i + 1, nm.c_str(), altloc, resnm.c_str(), ch, resnr,
 383                                      resic, 10 * x[i][XX], 10 * x[i][YY], 10 * x[i][ZZ], occup,
 384                                      bfac, atoms->atom[i].elem);
 385
 386             if (atoms->pdbinfo && atoms->pdbinfo[i].bAnisotropic)
 387             {
 388                 fprintf(out, "ANISOU%5d  %-4.4s%4.4s%c%4d%c %7d%7d%7d%7d%7d%7d\n", (i + 1) % 100000,
 389                         nm.c_str(), resnm.c_str(), ch, resnr, (resic == '\0') ? ' ' : resic,
 390                         atoms->pdbinfo[i].uij[0], atoms->pdbinfo[i].uij[1], atoms->pdbinfo[i].uij[2],
 391                         atoms->pdbinfo[i].uij[3], atoms->pdbinfo[i].uij[4], atoms->pdbinfo[i].uij[5]);
 392             }
 393         }
 394         else
 395         {
 396             gmx_fprintf_pqr_atomline(out, type, i + 1, nm.c_str(), resnm.c_str(), ch, resnr,
 397                                      10 * x[i][XX], 10 * x[i][YY], 10 * x[i][ZZ], occup, bfac);
 398         }
 399     }
 400
 401     fprintf(out, "TER\n");
 402     fprintf(out, "ENDMDL\n");
 403
 404     if (nullptr != gc)
 405     {
 406         /* Write conect records */
 407         for (int i = 0; (i < gc->nconect); i++)
 408         {
 409             fprintf(out, "CONECT%5d%5d\n", gc->conect[i].ai + 1, gc->conect[i].aj + 1);
 410         }
 411     }
 412 }
 413
 414 void write_pdbfile(FILE*          out,
 415                    const char*    title,
 416                    const t_atoms* atoms,
 417                    const rvec     x[],
 418                    int            ePBC,
 419                    const matrix   box,
 420                    char           chainid,
 421                    int            model_nr,
 422                    gmx_conect     conect)
 423 {
 424     int i, *index;
 425
 426     snew(index, atoms->nr);
 427     for (i = 0; i < atoms->nr; i++)
 428     {
 429         index[i] = i;
 430     }
 431     write_pdbfile_indexed(out, title, atoms, x, ePBC, box, chainid, model_nr, atoms->nr, index,
 432                           conect, false);
 433     sfree(index);
 434 }
 435
 436 static int line2type(const char* line)
 437 {
 438     int  k;
 439     char type[8];
 440
 441     for (k = 0; (k < 6); k++)
 442     {
 443         type[k] = line[k];
 444     }
 445     type[k] = '\0';
 446
 447     for (k = 0; (k < epdbNR); k++)
 448     {
 449         if (std::strncmp(type, pdbtp[k], strlen(pdbtp[k])) == 0)
 450         {
 451             break;
 452         }
 453     }
 454
 455     return k;
 456 }
 457
 458 static void read_anisou(char line[], int natom, t_atoms* atoms)
 459 {
 460     int  i, j, k, atomnr;
 461     char nc = '\0';
 462     char anr[12], anm[12];
 463
 464     /* Skip over type */
 465     j = 6;
 466     for (k = 0; (k < 5); k++, j++)
 467     {
 468         anr[k] = line[j];
 469     }
 470     anr[k] = nc;
 471     j++;
 472     for (k = 0; (k < 4); k++, j++)
 473     {
 474         anm[k] = line[j];
 475     }
 476     anm[k] = nc;
 477     j++;
 478
 479     /* Strip off spaces */
 480     trim(anm);
 481
 482     /* Search backwards for number and name only */
 483     atomnr = std::strtol(anr, nullptr, 10);
 484     for (i = natom - 1; (i >= 0); i--)
 485     {
 486         if ((std::strcmp(anm, *(atoms->atomname[i])) == 0) && (atomnr == atoms->pdbinfo[i].atomnr))
 487         {
 488             break;
 489         }
 490     }
 491     if (i < 0)
 492     {
 493         fprintf(stderr, "Skipping ANISOU record (atom %s %d not found)\n", anm, atomnr);
 494     }
 495     else
 496     {
 497         if (sscanf(line + 29, "%d%d%d%d%d%d", &atoms->pdbinfo[i].uij[U11], &atoms->pdbinfo[i].uij[U22],
 498                    &atoms->pdbinfo[i].uij[U33], &atoms->pdbinfo[i].uij[U12],
 499                    &atoms->pdbinfo[i].uij[U13], &atoms->pdbinfo[i].uij[U23])
 500             == 6)
 501         {
 502             atoms->pdbinfo[i].bAnisotropic = TRUE;
 503         }
 504         else
 505         {
 506             fprintf(stderr, "Invalid ANISOU record for atom %d\n", i);
 507             atoms->pdbinfo[i].bAnisotropic = FALSE;
 508         }
 509     }
 510 }
 511
 512 void get_pdb_atomnumber(const t_atoms* atoms, AtomProperties* aps)
 513 {
 514     int    i, atomnumber, len;
 515     size_t k;
 516     char   anm[6], anm_copy[6];
 517     char   nc = '\0';
 518     real   eval;
 519
 520     if (!atoms->pdbinfo)
 521     {
 522         gmx_incons("Trying to deduce atomnumbers when no pdb information is present");
 523     }
 524     for (i = 0; (i < atoms->nr); i++)
 525     {
 526         std::strcpy(anm, atoms->pdbinfo[i].atomnm);
 527         std::strcpy(anm_copy, atoms->pdbinfo[i].atomnm);
 528         bool atomNumberSet = false;
 529         len                = strlen(anm);
 530         if ((anm[0] != ' ') && ((len <= 2) || !std::isdigit(anm[2])))
 531         {
 532             anm_copy[2] = nc;
 533             if (aps->setAtomProperty(epropElement, "???", anm_copy, &eval))
 534             {
 535                 atomnumber    = gmx::roundToInt(eval);
 536                 atomNumberSet = true;
 537             }
 538             else
 539             {
 540                 anm_copy[1] = nc;
 541                 if (aps->setAtomProperty(epropElement, "???", anm_copy, &eval))
 542                 {
 543                     atomnumber    = gmx::roundToInt(eval);
 544                     atomNumberSet = true;
 545                 }
 546             }
 547         }
 548         if (!atomNumberSet)
 549         {
 550             k = 0;
 551             while ((k < std::strlen(anm)) && (std::isspace(anm[k]) || std::isdigit(anm[k])))
 552             {
 553                 k++;
 554             }
 555             anm_copy[0] = anm[k];
 556             anm_copy[1] = nc;
 557             if (aps->setAtomProperty(epropElement, "???", anm_copy, &eval))
 558             {
 559                 atomnumber    = gmx::roundToInt(eval);
 560                 atomNumberSet = true;
 561             }
 562         }
 563         std::string buf;
 564         if (atomNumberSet)
 565         {
 566             atoms->atom[i].atomnumber = atomnumber;
 567             buf                       = aps->elementFromAtomNumber(atomnumber);
 568             if (debug)
 569             {
 570                 fprintf(debug, "Atomnumber for atom '%s' is %d\n", anm, atomnumber);
 571             }
 572         }
 573         buf.resize(3);
 574         std::strncpy(atoms->atom[i].elem, buf.c_str(), 4);
 575     }
 576 }
 577
 578 static int
 579 read_atom(t_symtab* symtab, const char line[], int type, int natom, 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)", natom + 1, atoms->nr);
 593     }
 594
 595     /* Skip over type */
 596     j = 6;
 597     for (k = 0; (k < 5); k++, j++)
 598     {
 599         anr[k] = line[j];
 600     }
 601     anr[k] = nc;
 602     trim(anr);
 603     j++;
 604     for (k = 0; (k < 4); k++, j++)
 605     {
 606         anm[k] = line[j];
 607     }
 608     anm[k] = nc;
 609     std::strcpy(anm_copy, anm);
 610     rtrim(anm_copy);
 611     atomnumber = 0;
 612     trim(anm);
 613     altloc = line[j];
 614     j++;
 615     for (k = 0; (k < 4); k++, j++)
 616     {
 617         resnm[k] = line[j];
 618     }
 619     resnm[k] = nc;
 620     trim(resnm);
 621
 622     chainid = line[j];
 623     j++;
 624
 625     for (k = 0; (k < 4); k++, j++)
 626     {
 627         rnr[k] = line[j];
 628     }
 629     rnr[k] = nc;
 630     trim(rnr);
 631     resnr = std::strtol(rnr, nullptr, 10);
 632     resic = line[j];
 633     j += 4;
 634
 635     /* X,Y,Z Coordinate */
 636     for (k = 0; (k < 8); k++, j++)
 637     {
 638         xc[k] = line[j];
 639     }
 640     xc[k] = nc;
 641     for (k = 0; (k < 8); k++, j++)
 642     {
 643         yc[k] = line[j];
 644     }
 645     yc[k] = nc;
 646     for (k = 0; (k < 8); k++, j++)
 647     {
 648         zc[k] = line[j];
 649     }
 650     zc[k] = nc;
 651
 652     /* Weight */
 653     for (k = 0; (k < 6); k++, j++)
 654     {
 655         occup[k] = line[j];
 656     }
 657     occup[k] = nc;
 658
 659     /* B-Factor */
 660     for (k = 0; (k < 7); k++, j++)
 661     {
 662         bfac[k] = line[j];
 663     }
 664     bfac[k] = nc;
 665
 666     /* 10 blanks */
 667     j += 10;
 668
 669     /* Element name */
 670     for (k = 0; (k < 2); k++, j++)
 671     {
 672         elem[k] = line[j];
 673     }
 674     elem[k] = nc;
 675     trim(elem);
 676
 677     if (atoms->atom)
 678     {
 679         atomn = &(atoms->atom[natom]);
 680         if ((natom == 0) || atoms->resinfo[atoms->atom[natom - 1].resind].nr != resnr
 681             || atoms->resinfo[atoms->atom[natom - 1].resind].ic != resic
 682             || (strcmp(*atoms->resinfo[atoms->atom[natom - 1].resind].name, resnm) != 0))
 683         {
 684             if (natom == 0)
 685             {
 686                 atomn->resind = 0;
 687             }
 688             else
 689             {
 690                 atomn->resind = atoms->atom[natom - 1].resind + 1;
 691             }
 692             atoms->nres = atomn->resind + 1;
 693             t_atoms_set_resinfo(atoms, natom, symtab, resnm, resnr, resic, chainnum, chainid);
 694         }
 695         else
 696         {
 697             atomn->resind = atoms->atom[natom - 1].resind;
 698         }
 699         if (bChange)
 700         {
 701             xlate_atomname_pdb2gmx(anm);
 702         }
 703         atoms->atomname[natom] = put_symtab(symtab, anm);
 704         atomn->m               = 0.0;
 705         atomn->q               = 0.0;
 706         atomn->atomnumber      = atomnumber;
 707         strncpy(atomn->elem, elem, 4);
 708     }
 709     x[natom][XX] = strtod(xc, nullptr) * 0.1;
 710     x[natom][YY] = strtod(yc, nullptr) * 0.1;
 711     x[natom][ZZ] = strtod(zc, nullptr) * 0.1;
 712     if (atoms->pdbinfo)
 713     {
 714         atoms->pdbinfo[natom].type   = type;
 715         atoms->pdbinfo[natom].atomnr = strtol(anr, nullptr, 10);
 716         atoms->pdbinfo[natom].altloc = altloc;
 717         strcpy(atoms->pdbinfo[natom].atomnm, anm_copy);
 718         atoms->pdbinfo[natom].bfac  = strtod(bfac, nullptr);
 719         atoms->pdbinfo[natom].occup = strtod(occup, nullptr);
 720     }
 721     natom++;
 722
 723     return natom;
 724 }
 725
 726 gmx_bool is_hydrogen(const char* nm)
 727 {
 728     char buf[30];
 729
 730     std::strcpy(buf, nm);
 731     trim(buf);
 732
 733     if (buf[0] == 'H')
 734     {
 735         return TRUE;
 736     }
 737     else if ((std::isdigit(buf[0])) && (buf[1] == 'H'))
 738     {
 739         return TRUE;
 740     }
 741     return FALSE;
 742 }
 743
 744 gmx_bool is_dummymass(const char* nm)
 745 {
 746     char buf[30];
 747
 748     std::strcpy(buf, nm);
 749     trim(buf);
 750
 751     return (buf[0] == 'M') && (std::isdigit(buf[strlen(buf) - 1]) != 0);
 752 }
 753
 754 static void gmx_conect_addline(gmx_conect_t* con, char* line)
 755 {
 756     int n, ai, aj;
 757
 758     std::string form2  = "%*s";
 759     std::string format = form2 + "%d";
 760     if (sscanf(line, format.c_str(), &ai) == 1)
 761     {
 762         do
 763         {
 764             form2 += "%*s";
 765             format = form2 + "%d";
 766             n      = sscanf(line, format.c_str(), &aj);
 767             if (n == 1)
 768             {
 769                 gmx_conect_add(con, ai - 1, aj - 1); /* to prevent duplicated records */
 770             }
 771         } while (n == 1);
 772     }
 773 }
 774
 775 void gmx_conect_dump(FILE* fp, gmx_conect conect)
 776 {
 777     gmx_conect_t* gc = static_cast<gmx_conect_t*>(conect);
 778     int           i;
 779
 780     for (i = 0; (i < gc->nconect); i++)
 781     {
 782         fprintf(fp, "%6s%5d%5d\n", "CONECT", gc->conect[i].ai + 1, gc->conect[i].aj + 1);
 783     }
 784 }
 785
 786 gmx_conect gmx_conect_init()
 787 {
 788     gmx_conect_t* gc;
 789
 790     snew(gc, 1);
 791
 792     return gc;
 793 }
 794
 795 void gmx_conect_done(gmx_conect conect)
 796 {
 797     gmx_conect_t* gc = conect;
 798
 799     sfree(gc->conect);
 800 }
 801
 802 gmx_bool gmx_conect_exist(gmx_conect conect, int ai, int aj)
 803 {
 804     gmx_conect_t* gc = conect;
 805     int           i;
 806
 807     /* if (!gc->bSorted)
 808        sort_conect(gc);*/
 809
 810     for (i = 0; (i < gc->nconect); i++)
 811     {
 812         if (((gc->conect[i].ai == ai) && (gc->conect[i].aj == aj))
 813             || ((gc->conect[i].aj == ai) && (gc->conect[i].ai == aj)))
 814         {
 815             return TRUE;
 816         }
 817     }
 818     return FALSE;
 819 }
 820
 821 void gmx_conect_add(gmx_conect conect, int ai, int aj)
 822 {
 823     gmx_conect_t* gc = static_cast<gmx_conect_t*>(conect);
 824
 825     /* if (!gc->bSorted)
 826        sort_conect(gc);*/
 827
 828     if (!gmx_conect_exist(conect, ai, aj))
 829     {
 830         srenew(gc->conect, ++gc->nconect);
 831         gc->conect[gc->nconect - 1].ai = ai;
 832         gc->conect[gc->nconect - 1].aj = aj;
 833     }
 834 }
 835
 836 int read_pdbfile(FILE*      in,
 837                  char*      title,
 838                  int*       model_nr,
 839                  t_atoms*   atoms,
 840                  t_symtab*  symtab,
 841                  rvec       x[],
 842                  int*       ePBC,
 843                  matrix     box,
 844                  gmx_bool   bChange,
 845                  gmx_conect conect)
 846 {
 847     gmx_conect_t* gc = conect;
 848     gmx_bool      bCOMPND;
 849     gmx_bool      bConnWarn = FALSE;
 850     char          line[STRLEN + 1];
 851     int           line_type;
 852     char *        c, *d;
 853     int           natom, chainnum;
 854     gmx_bool      bStop = FALSE;
 855
 856     if (ePBC)
 857     {
 858         /* Only assume pbc when there is a CRYST1 entry */
 859         *ePBC = epbcNONE;
 860     }
 861     if (box != nullptr)
 862     {
 863         clear_mat(box);
 864     }
 865
 866     atoms->haveMass    = FALSE;
 867     atoms->haveCharge  = FALSE;
 868     atoms->haveType    = FALSE;
 869     atoms->haveBState  = FALSE;
 870     atoms->havePdbInfo = (atoms->pdbinfo != nullptr);
 871
 872     bCOMPND  = FALSE;
 873     title[0] = '\0';
 874     natom    = 0;
 875     chainnum = 0;
 876     while (!bStop && (fgets2(line, STRLEN, in) != nullptr))
 877     {
 878         line_type = line2type(line);
 879
 880         switch (line_type)
 881         {
 882             case epdbATOM:
 883             case epdbHETATM:
 884                 natom = read_atom(symtab, line, line_type, natom, atoms, x, chainnum, bChange);
 885                 break;
 886
 887             case epdbANISOU:
 888                 if (atoms->havePdbInfo)
 889                 {
 890                     read_anisou(line, natom, atoms);
 891                 }
 892                 break;
 893
 894             case epdbCRYST1: read_cryst1(line, ePBC, box); break;
 895
 896             case epdbTITLE:
 897             case epdbHEADER:
 898                 if (std::strlen(line) > 6)
 899                 {
 900                     c = line + 6;
 901                     /* skip HEADER or TITLE and spaces */
 902                     while (c[0] != ' ')
 903                     {
 904                         c++;
 905                     }
 906                     while (c[0] == ' ')
 907                     {
 908                         c++;
 909                     }
 910                     /* truncate after title */
 911                     d = std::strstr(c, "      ");
 912                     if (d)
 913                     {
 914                         d[0] = '\0';
 915                     }
 916                     if (std::strlen(c) > 0)
 917                     {
 918                         std::strcpy(title, c);
 919                     }
 920                 }
 921                 break;
 922
 923             case epdbCOMPND:
 924                 if ((!std::strstr(line, ": ")) || (std::strstr(line + 6, "MOLECULE:")))
 925                 {
 926                     if (!(c = std::strstr(line + 6, "MOLECULE:")))
 927                     {
 928                         c = line;
 929                     }
 930                     /* skip 'MOLECULE:' and spaces */
 931                     while (c[0] != ' ')
 932                     {
 933                         c++;
 934                     }
 935                     while (c[0] == ' ')
 936                     {
 937                         c++;
 938                     }
 939                     /* truncate after title */
 940                     d = strstr(c, "   ");
 941                     if (d)
 942                     {
 943                         while ((d[-1] == ';') && d > c)
 944                         {
 945                             d--;
 946                         }
 947                         d[0] = '\0';
 948                     }
 949                     if (strlen(c) > 0)
 950                     {
 951                         if (bCOMPND)
 952                         {
 953                             std::strcat(title, "; ");
 954                             std::strcat(title, c);
 955                         }
 956                         else
 957                         {
 958                             std::strcpy(title, c);
 959                         }
 960                     }
 961                     bCOMPND = TRUE;
 962                 }
 963                 break;
 964
 965             case epdbTER: chainnum++; break;
 966
 967             case epdbMODEL:
 968                 if (model_nr)
 969                 {
 970                     sscanf(line, "%*s%d", model_nr);
 971                 }
 972                 break;
 973
 974             case epdbENDMDL: bStop = TRUE; break;
 975             case epdbCONECT:
 976                 if (gc)
 977                 {
 978                     gmx_conect_addline(gc, line);
 979                 }
 980                 else if (!bConnWarn)
 981                 {
 982                     fprintf(stderr, "WARNING: all CONECT records are ignored\n");
 983                     bConnWarn = TRUE;
 984                 }
 985                 break;
 986
 987             default: break;
 988         }
 989     }
 990
 991     return natom;
 992 }
 993
 994 void get_pdb_coordnum(FILE* in, int* natoms)
 995 {
 996     char line[STRLEN];
 997
 998     *natoms = 0;
 999     while (fgets2(line, STRLEN, in))
1000     {
1001         if (std::strncmp(line, "ENDMDL", 6) == 0)
1002         {
1003             break;
1004         }
1005         if ((std::strncmp(line, "ATOM  ", 6) == 0) || (std::strncmp(line, "HETATM", 6) == 0))
1006         {
1007             (*natoms)++;
1008         }
1009     }
1010 }
1011
1012 void gmx_pdb_read_conf(const char* infile, t_symtab* symtab, char** name, t_atoms* atoms, rvec x[], int* ePBC, matrix box)
1013 {
1014     FILE* in = gmx_fio_fopen(infile, "r");
1015     char  title[STRLEN];
1016     read_pdbfile(in, title, nullptr, atoms, symtab, x, ePBC, box, TRUE, nullptr);
1017     if (name != nullptr)
1018     {
1019         *name = gmx_strdup(title);
1020     }
1021     gmx_fio_fclose(in);
1022 }
1023
1024 gmx_conect gmx_conect_generate(const t_topology* top)
1025 {
1026     int        f, i;
1027     gmx_conect gc;
1028
1029     /* Fill the conect records */
1030     gc = gmx_conect_init();
1031
1032     for (f = 0; (f < F_NRE); f++)
1033     {
1034         if (IS_CHEMBOND(f))
1035         {
1036             for (i = 0; (i < top->idef.il[f].nr); i += interaction_function[f].nratoms + 1)
1037             {
1038                 gmx_conect_add(gc, top->idef.il[f].iatoms[i + 1], top->idef.il[f].iatoms[i + 2]);
1039             }
1040         }
1041     }
1042     return gc;
1043 }
1044
1045 int gmx_fprintf_pdb_atomline(FILE*           fp,
1046                              enum PDB_record record,
1047                              int             atom_seq_number,
1048                              const char*     atom_name,
1049                              char            alternate_location,
1050                              const char*     res_name,
1051                              char            chain_id,
1052                              int             res_seq_number,
1053                              char            res_insertion_code,
1054                              real            x,
1055                              real            y,
1056                              real            z,
1057                              real            occupancy,
1058                              real            b_factor,
1059                              const char*     element)
1060 {
1061     char     tmp_atomname[6], tmp_resname[6];
1062     gmx_bool start_name_in_col13;
1063     int      n;
1064
1065     if (record != epdbATOM && record != epdbHETATM)
1066     {
1067         gmx_fatal(FARGS, "Can only print PDB atom lines as ATOM or HETATM records");
1068     }
1069
1070     /* Format atom name */
1071     if (atom_name != nullptr)
1072     {
1073         /* If the atom name is an element name with two chars, it should start already in column 13.
1074          * Otherwise it should start in column 14, unless the name length is 4 chars.
1075          */
1076         if ((element != nullptr) && (std::strlen(element) >= 2)
1077             && (gmx_strncasecmp(atom_name, element, 2) == 0))
1078         {
1079             start_name_in_col13 = TRUE;
1080         }
1081         else
1082         {
1083             start_name_in_col13 = (std::strlen(atom_name) >= 4);
1084         }
1085         snprintf(tmp_atomname, sizeof(tmp_atomname), start_name_in_col13 ? "" : " ");
1086         std::strncat(tmp_atomname, atom_name, 4);
1087         tmp_atomname[5] = '\0';
1088     }
1089     else
1090     {
1091         tmp_atomname[0] = '\0';
1092     }
1093
1094     /* Format residue name */
1095     std::strncpy(tmp_resname, (res_name != nullptr) ? res_name : "", 4);
1096     /* Make sure the string is terminated if strlen was > 4 */
1097     tmp_resname[4] = '\0';
1098     /* String is properly terminated, so now we can use strcat. By adding a
1099      * space we can write it right-justified, and if the original name was
1100      * three characters or less there will be a space added on the right side.
1101      */
1102     std::strcat(tmp_resname, " ");
1103
1104     /* Truncate integers so they fit */
1105     atom_seq_number = atom_seq_number % 100000;
1106     res_seq_number  = res_seq_number % 10000;
1107
1108     n = fprintf(fp, "%-6s%5d %-4.4s%c%4.4s%c%4d%c   %8.3f%8.3f%8.3f%6.2f%6.2f          %2s\n",
1109                 pdbtp[record], atom_seq_number, tmp_atomname, alternate_location, tmp_resname,
1110                 chain_id, res_seq_number, res_insertion_code, x, y, z, occupancy, b_factor,
1111                 (element != nullptr) ? element : "");
1112
1113     return n;
1114 }