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,2021, 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 "gromacs/fileio/pdbio.h"
  41
  42 #include <algorithm>
  43 #include <cctype>
  44 #include <cmath>
  45 #include <cstdlib>
  46 #include <cstring>
  47
  48 #include <string>
  49
  50 #include "gromacs/fileio/gmxfio.h"
  51 #include "gromacs/math/units.h"
  52 #include "gromacs/math/utilities.h"
  53 #include "gromacs/math/vec.h"
  54 #include "gromacs/pbcutil/pbc.h"
  55 #include "gromacs/topology/atomprop.h"
  56 #include "gromacs/topology/ifunc.h"
  57 #include "gromacs/topology/residuetypes.h"
  58 #include "gromacs/topology/symtab.h"
  59 #include "gromacs/topology/topology.h"
  60 #include "gromacs/utility/coolstuff.h"
  61 #include "gromacs/utility/cstringutil.h"
  62 #include "gromacs/utility/enumerationhelpers.h"
  63 #include "gromacs/utility/fatalerror.h"
  64 #include "gromacs/utility/futil.h"
  65 #include "gromacs/utility/smalloc.h"
  66 #include "gromacs/utility/snprintf.h"
  67 #include "gromacs/utility/stringtoenumvalueconverter.h"
  68
  69 typedef struct
  70 {
  71     int ai, aj;
  72 } gmx_conection_t;
  73
  74 typedef struct gmx_conect_t
  75 {
  76     int              nconect;
  77     gmx_bool         bSorted;
  78     gmx_conection_t* conect;
  79 } gmx_conect_t;
  80
  81 const char* enumValueToString(PdbRecordType enumValue)
  82 {
  83     static constexpr gmx::EnumerationArray<PdbRecordType, const char*> pdbRecordTypeName = {
  84         "ATOM  ", "HETATM", "ANISOU", "CRYST1", "COMPND", "MODEL",
  85         "ENDMDL", "TER",    "HEADER", "TITLE",  "REMARK", "CONECT"
  86     };
  87     return pdbRecordTypeName[enumValue];
  88 }
  89
  90 void gmx_write_pdb_box(FILE* out, PbcType pbcType, const matrix box)
  91 {
  92     real alpha, beta, gamma;
  93
  94     if (pbcType == PbcType::Unset)
  95     {
  96         pbcType = guessPbcType(box);
  97     }
  98
  99     if (pbcType == PbcType::No)
 100     {
 101         return;
 102     }
 103
 104     if (norm2(box[YY]) * norm2(box[ZZ]) != 0)
 105     {
 106         alpha = gmx::c_rad2Deg * gmx_angle(box[YY], box[ZZ]);
 107     }
 108     else
 109     {
 110         alpha = 90;
 111     }
 112     if (norm2(box[XX]) * norm2(box[ZZ]) != 0)
 113     {
 114         beta = gmx::c_rad2Deg * gmx_angle(box[XX], box[ZZ]);
 115     }
 116     else
 117     {
 118         beta = 90;
 119     }
 120     if (norm2(box[XX]) * norm2(box[YY]) != 0)
 121     {
 122         gamma = gmx::c_rad2Deg * gmx_angle(box[XX], box[YY]);
 123     }
 124     else
 125     {
 126         gamma = 90;
 127     }
 128     fprintf(out, "REMARK    THIS IS A SIMULATION BOX\n");
 129     if (pbcType != PbcType::Screw)
 130     {
 131         fprintf(out,
 132                 "CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f %-11s%4d\n",
 133                 10 * norm(box[XX]),
 134                 10 * norm(box[YY]),
 135                 10 * norm(box[ZZ]),
 136                 alpha,
 137                 beta,
 138                 gamma,
 139                 "P 1",
 140                 1);
 141     }
 142     else
 143     {
 144         /* Double the a-vector length and write the correct space group */
 145         fprintf(out,
 146                 "CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f %-11s%4d\n",
 147                 20 * norm(box[XX]),
 148                 10 * norm(box[YY]),
 149                 10 * norm(box[ZZ]),
 150                 alpha,
 151                 beta,
 152                 gamma,
 153                 "P 21 1 1",
 154                 1);
 155     }
 156 }
 157
 158 static void read_cryst1(char* line, PbcType* pbcType, matrix box)
 159 {
 160 #define SG_SIZE 11
 161     char    sa[12], sb[12], sc[12], sg[SG_SIZE + 1], ident;
 162     double  fa, fb, fc, alpha, beta, gamma, cosa, cosb, cosg, sing;
 163     int     syma, symb, symc;
 164     PbcType pbcTypeFile;
 165
 166     sscanf(line, "%*s%s%s%s%lf%lf%lf", sa, sb, sc, &alpha, &beta, &gamma);
 167
 168     pbcTypeFile = PbcType::Unset;
 169     if (strlen(line) >= 55)
 170     {
 171         strncpy(sg, line + 55, SG_SIZE);
 172         sg[SG_SIZE] = '\0';
 173         ident       = ' ';
 174         syma        = 0;
 175         symb        = 0;
 176         symc        = 0;
 177         sscanf(sg, "%c %d %d %d", &ident, &syma, &symb, &symc);
 178         if (ident == 'P' && syma == 1 && symb <= 1 && symc <= 1)
 179         {
 180             fc          = strtod(sc, nullptr) * 0.1;
 181             pbcTypeFile = (fc > 0 ? PbcType::Xyz : PbcType::XY);
 182         }
 183         if (ident == 'P' && syma == 21 && symb == 1 && symc == 1)
 184         {
 185             pbcTypeFile = PbcType::Screw;
 186         }
 187     }
 188     if (pbcType)
 189     {
 190         *pbcType = pbcTypeFile;
 191     }
 192
 193     if (box)
 194     {
 195         fa = strtod(sa, nullptr) * 0.1;
 196         fb = strtod(sb, nullptr) * 0.1;
 197         fc = strtod(sc, nullptr) * 0.1;
 198         if (pbcTypeFile == PbcType::Screw)
 199         {
 200             fa *= 0.5;
 201         }
 202         clear_mat(box);
 203         box[XX][XX] = fa;
 204         if ((alpha != 90.0) || (beta != 90.0) || (gamma != 90.0))
 205         {
 206             if (alpha != 90.0)
 207             {
 208                 cosa = std::cos(alpha * gmx::c_deg2Rad);
 209             }
 210             else
 211             {
 212                 cosa = 0;
 213             }
 214             if (beta != 90.0)
 215             {
 216                 cosb = std::cos(beta * gmx::c_deg2Rad);
 217             }
 218             else
 219             {
 220                 cosb = 0;
 221             }
 222             if (gamma != 90.0)
 223             {
 224                 cosg = std::cos(gamma * gmx::c_deg2Rad);
 225                 sing = std::sin(gamma * gmx::c_deg2Rad);
 226             }
 227             else
 228             {
 229                 cosg = 0;
 230                 sing = 1;
 231             }
 232             box[YY][XX] = fb * cosg;
 233             box[YY][YY] = fb * sing;
 234             box[ZZ][XX] = fc * cosb;
 235             box[ZZ][YY] = fc * (cosa - cosb * cosg) / sing;
 236             box[ZZ][ZZ] = std::sqrt(fc * fc - box[ZZ][XX] * box[ZZ][XX] - box[ZZ][YY] * box[ZZ][YY]);
 237         }
 238         else
 239         {
 240             box[YY][YY] = fb;
 241             box[ZZ][ZZ] = fc;
 242         }
 243     }
 244 }
 245
 246 static int gmx_fprintf_pqr_atomline(FILE*         fp,
 247                                     PdbRecordType record,
 248                                     int           atom_seq_number,
 249                                     const char*   atom_name,
 250                                     const char*   res_name,
 251                                     char          chain_id,
 252                                     int           res_seq_number,
 253                                     real          x,
 254                                     real          y,
 255                                     real          z,
 256                                     real          occupancy,
 257                                     real          b_factor)
 258 {
 259     GMX_RELEASE_ASSERT(record == PdbRecordType::Atom || record == PdbRecordType::Hetatm,
 260                        "Can only print PQR atom lines as ATOM or HETATM records");
 261
 262     /* Check atom name */
 263     GMX_RELEASE_ASSERT(atom_name != nullptr, "Need atom information to print pqr");
 264
 265     /* Check residue name */
 266     GMX_RELEASE_ASSERT(res_name != nullptr, "Need residue information to print pqr");
 267
 268     /* Truncate integers so they fit */
 269     atom_seq_number = atom_seq_number % 100000;
 270     res_seq_number  = res_seq_number % 10000;
 271
 272     int n = fprintf(fp,
 273                     "%-6s%5d %-4.4s%4.4s%c%4d %8.3f %8.3f %8.3f %6.2f %6.2f\n",
 274                     enumValueToString(record),
 275                     atom_seq_number,
 276                     atom_name,
 277                     res_name,
 278                     chain_id,
 279                     res_seq_number,
 280                     x,
 281                     y,
 282                     z,
 283                     occupancy,
 284                     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                            PbcType        pbcType,
 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     PdbRecordType 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, pbcType, 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         int           resnr = atoms->resinfo[resind].nr;
 342         unsigned char resic = atoms->resinfo[resind].ic;
 343         unsigned char ch;
 344         if (chainid != ' ')
 345         {
 346             ch = chainid;
 347         }
 348         else
 349         {
 350             ch = atoms->resinfo[resind].chainid;
 351
 352             if (ch == 0)
 353             {
 354                 ch = ' ';
 355             }
 356         }
 357         if (resnr >= 10000)
 358         {
 359             resnr = resnr % 10000;
 360         }
 361         t_pdbinfo pdbinfo;
 362         if (atoms->pdbinfo != nullptr)
 363         {
 364             pdbinfo = atoms->pdbinfo[i];
 365         }
 366         else
 367         {
 368             gmx_pdbinfo_init_default(&pdbinfo);
 369         }
 370         type   = pdbinfo.type;
 371         altloc = pdbinfo.altloc;
 372         if (!isalnum(altloc))
 373         {
 374             altloc = ' ';
 375         }
 376         occup = bOccup ? 1.0 : pdbinfo.occup;
 377         bfac  = pdbinfo.bfac;
 378         if (!usePqrFormat)
 379         {
 380             gmx_fprintf_pdb_atomline(out,
 381                                      type,
 382                                      i + 1,
 383                                      nm.c_str(),
 384                                      altloc,
 385                                      resnm.c_str(),
 386                                      ch,
 387                                      resnr,
 388                                      resic,
 389                                      10 * x[i][XX],
 390                                      10 * x[i][YY],
 391                                      10 * x[i][ZZ],
 392                                      occup,
 393                                      bfac,
 394                                      atoms->atom[i].elem);
 395
 396             if (atoms->pdbinfo && atoms->pdbinfo[i].bAnisotropic)
 397             {
 398                 fprintf(out,
 399                         "ANISOU%5d  %-4.4s%4.4s%c%4d%c %7d%7d%7d%7d%7d%7d\n",
 400                         (i + 1) % 100000,
 401                         nm.c_str(),
 402                         resnm.c_str(),
 403                         ch,
 404                         resnr,
 405                         (resic == '\0') ? ' ' : resic,
 406                         atoms->pdbinfo[i].uij[0],
 407                         atoms->pdbinfo[i].uij[1],
 408                         atoms->pdbinfo[i].uij[2],
 409                         atoms->pdbinfo[i].uij[3],
 410                         atoms->pdbinfo[i].uij[4],
 411                         atoms->pdbinfo[i].uij[5]);
 412             }
 413         }
 414         else
 415         {
 416             gmx_fprintf_pqr_atomline(out,
 417                                      type,
 418                                      i + 1,
 419                                      nm.c_str(),
 420                                      resnm.c_str(),
 421                                      ch,
 422                                      resnr,
 423                                      10 * x[i][XX],
 424                                      10 * x[i][YY],
 425                                      10 * x[i][ZZ],
 426                                      occup,
 427                                      bfac);
 428         }
 429     }
 430
 431     fprintf(out, "TER\n");
 432     fprintf(out, "ENDMDL\n");
 433
 434     if (nullptr != gc)
 435     {
 436         /* Write conect records */
 437         for (int i = 0; (i < gc->nconect); i++)
 438         {
 439             fprintf(out, "CONECT%5d%5d\n", gc->conect[i].ai + 1, gc->conect[i].aj + 1);
 440         }
 441     }
 442 }
 443
 444 void write_pdbfile(FILE*          out,
 445                    const char*    title,
 446                    const t_atoms* atoms,
 447                    const rvec     x[],
 448                    PbcType        pbcType,
 449                    const matrix   box,
 450                    char           chainid,
 451                    int            model_nr,
 452                    gmx_conect     conect)
 453 {
 454     int i, *index;
 455
 456     snew(index, atoms->nr);
 457     for (i = 0; i < atoms->nr; i++)
 458     {
 459         index[i] = i;
 460     }
 461     write_pdbfile_indexed(
 462             out, title, atoms, x, pbcType, box, chainid, model_nr, atoms->nr, index, conect, false);
 463     sfree(index);
 464 }
 465
 466 static void read_anisou(char line[], int natom, t_atoms* atoms)
 467 {
 468     int  i, j, k, atomnr;
 469     char nc = '\0';
 470     char anr[12], anm[12];
 471
 472     /* Skip over type */
 473     j = 6;
 474     for (k = 0; (k < 5); k++, j++)
 475     {
 476         anr[k] = line[j];
 477     }
 478     anr[k] = nc;
 479     j++;
 480     for (k = 0; (k < 4); k++, j++)
 481     {
 482         anm[k] = line[j];
 483     }
 484     anm[k] = nc;
 485     j++;
 486
 487     /* Strip off spaces */
 488     trim(anm);
 489
 490     /* Search backwards for number and name only */
 491     atomnr = std::strtol(anr, nullptr, 10);
 492     for (i = natom - 1; (i >= 0); i--)
 493     {
 494         if ((std::strcmp(anm, *(atoms->atomname[i])) == 0) && (atomnr == atoms->pdbinfo[i].atomnr))
 495         {
 496             break;
 497         }
 498     }
 499     if (i < 0)
 500     {
 501         fprintf(stderr, "Skipping ANISOU record (atom %s %d not found)\n", anm, atomnr);
 502     }
 503     else
 504     {
 505         if (sscanf(line + 29,
 506                    "%d%d%d%d%d%d",
 507                    &atoms->pdbinfo[i].uij[U11],
 508                    &atoms->pdbinfo[i].uij[U22],
 509                    &atoms->pdbinfo[i].uij[U33],
 510                    &atoms->pdbinfo[i].uij[U12],
 511                    &atoms->pdbinfo[i].uij[U13],
 512                    &atoms->pdbinfo[i].uij[U23])
 513             == 6)
 514         {
 515             atoms->pdbinfo[i].bAnisotropic = TRUE;
 516         }
 517         else
 518         {
 519             fprintf(stderr, "Invalid ANISOU record for atom %d\n", i);
 520             atoms->pdbinfo[i].bAnisotropic = FALSE;
 521         }
 522     }
 523 }
 524
 525 void get_pdb_atomnumber(const t_atoms* atoms, AtomProperties* aps)
 526 {
 527     int    i, atomnumber, len;
 528     size_t k;
 529     char   anm[6], anm_copy[6];
 530     char   nc = '\0';
 531     real   eval;
 532
 533     if (!atoms->pdbinfo)
 534     {
 535         gmx_incons("Trying to deduce atomnumbers when no pdb information is present");
 536     }
 537     for (i = 0; (i < atoms->nr); i++)
 538     {
 539         std::strcpy(anm, atoms->pdbinfo[i].atomnm);
 540         std::strcpy(anm_copy, atoms->pdbinfo[i].atomnm);
 541         bool atomNumberSet = false;
 542         len                = strlen(anm);
 543         if ((anm[0] != ' ') && ((len <= 2) || !std::isdigit(anm[2])))
 544         {
 545             anm_copy[2] = nc;
 546             if (aps->setAtomProperty(epropElement, "???", anm_copy, &eval))
 547             {
 548                 atomnumber    = gmx::roundToInt(eval);
 549                 atomNumberSet = true;
 550             }
 551             else
 552             {
 553                 anm_copy[1] = nc;
 554                 if (aps->setAtomProperty(epropElement, "???", anm_copy, &eval))
 555                 {
 556                     atomnumber    = gmx::roundToInt(eval);
 557                     atomNumberSet = true;
 558                 }
 559             }
 560         }
 561         if (!atomNumberSet)
 562         {
 563             k = 0;
 564             while ((k < std::strlen(anm)) && (std::isspace(anm[k]) || std::isdigit(anm[k])))
 565             {
 566                 k++;
 567             }
 568             anm_copy[0] = anm[k];
 569             anm_copy[1] = nc;
 570             if (aps->setAtomProperty(epropElement, "???", anm_copy, &eval))
 571             {
 572                 atomnumber    = gmx::roundToInt(eval);
 573                 atomNumberSet = true;
 574             }
 575         }
 576         static constexpr size_t sc_maxElementNameLength = 3;
 577         static_assert(sizeof(atoms->atom[i].elem) >= sc_maxElementNameLength + 1);
 578         std::string element;
 579         if (atomNumberSet)
 580         {
 581             atoms->atom[i].atomnumber = atomnumber;
 582             element                   = aps->elementFromAtomNumber(atomnumber);
 583             if (debug)
 584             {
 585                 fprintf(debug, "Atomnumber for atom '%s' is %d\n", anm, atomnumber);
 586             }
 587         }
 588         element.resize(sc_maxElementNameLength);
 589         std::strcpy(atoms->atom[i].elem, element.c_str());
 590     }
 591 }
 592
 593 static int
 594 read_atom(t_symtab* symtab, const char line[], PdbRecordType type, int natom, t_atoms* atoms, rvec x[], int chainnum)
 595 {
 596     t_atom*       atomn;
 597     int           j, k;
 598     char          nc = '\0';
 599     char          anr[12], anm[12], anm_copy[12], altloc, resnm[12], rnr[12], elem[3];
 600     char          xc[12], yc[12], zc[12], occup[12], bfac[12];
 601     unsigned char resic;
 602     char          chainid;
 603     int           resnr, atomnumber;
 604
 605     if (natom >= atoms->nr)
 606     {
 607         gmx_fatal(FARGS, "\nFound more atoms (%d) in pdb file than expected (%d)", natom + 1, atoms->nr);
 608     }
 609
 610     /* Skip over type */
 611     j = 6;
 612     for (k = 0; (k < 5); k++, j++)
 613     {
 614         anr[k] = line[j];
 615     }
 616     anr[k] = nc;
 617     trim(anr);
 618     j++;
 619     for (k = 0; (k < 4); k++, j++)
 620     {
 621         anm[k] = line[j];
 622     }
 623     anm[k] = nc;
 624     std::strcpy(anm_copy, anm);
 625     rtrim(anm_copy);
 626     atomnumber = 0;
 627     trim(anm);
 628     altloc = line[j];
 629     j++;
 630     for (k = 0; (k < 4); k++, j++)
 631     {
 632         resnm[k] = line[j];
 633     }
 634     resnm[k] = nc;
 635     trim(resnm);
 636
 637     chainid = line[j];
 638     j++;
 639
 640     for (k = 0; (k < 4); k++, j++)
 641     {
 642         rnr[k] = line[j];
 643     }
 644     rnr[k] = nc;
 645     trim(rnr);
 646     resnr = std::strtol(rnr, nullptr, 10);
 647     resic = line[j];
 648     j += 4;
 649
 650     /* X,Y,Z Coordinate */
 651     for (k = 0; (k < 8); k++, j++)
 652     {
 653         xc[k] = line[j];
 654     }
 655     xc[k] = nc;
 656     for (k = 0; (k < 8); k++, j++)
 657     {
 658         yc[k] = line[j];
 659     }
 660     yc[k] = nc;
 661     for (k = 0; (k < 8); k++, j++)
 662     {
 663         zc[k] = line[j];
 664     }
 665     zc[k] = nc;
 666
 667     /* Weight */
 668     for (k = 0; (k < 6); k++, j++)
 669     {
 670         occup[k] = line[j];
 671     }
 672     occup[k] = nc;
 673
 674     /* B-Factor */
 675     for (k = 0; (k < 7); k++, j++)
 676     {
 677         bfac[k] = line[j];
 678     }
 679     bfac[k] = nc;
 680
 681     /* 10 blanks */
 682     j += 10;
 683
 684     /* Element name */
 685     for (k = 0; (k < 2); k++, j++)
 686     {
 687         elem[k] = line[j];
 688     }
 689     elem[k] = nc;
 690     trim(elem);
 691
 692     if (atoms->atom)
 693     {
 694         atomn = &(atoms->atom[natom]);
 695         if ((natom == 0) || atoms->resinfo[atoms->atom[natom - 1].resind].nr != resnr
 696             || atoms->resinfo[atoms->atom[natom - 1].resind].ic != resic
 697             || (strcmp(*atoms->resinfo[atoms->atom[natom - 1].resind].name, resnm) != 0))
 698         {
 699             if (natom == 0)
 700             {
 701                 atomn->resind = 0;
 702             }
 703             else
 704             {
 705                 atomn->resind = atoms->atom[natom - 1].resind + 1;
 706             }
 707             atoms->nres = atomn->resind + 1;
 708             t_atoms_set_resinfo(atoms, natom, symtab, resnm, resnr, resic, chainnum, chainid);
 709         }
 710         else
 711         {
 712             atomn->resind = atoms->atom[natom - 1].resind;
 713         }
 714         atoms->atomname[natom] = put_symtab(symtab, anm);
 715         atomn->m               = 0.0;
 716         atomn->q               = 0.0;
 717         atomn->atomnumber      = atomnumber;
 718         strncpy(atomn->elem, elem, 4);
 719     }
 720     x[natom][XX] = strtod(xc, nullptr) * 0.1;
 721     x[natom][YY] = strtod(yc, nullptr) * 0.1;
 722     x[natom][ZZ] = strtod(zc, nullptr) * 0.1;
 723     if (atoms->pdbinfo)
 724     {
 725         atoms->pdbinfo[natom].type   = type;
 726         atoms->pdbinfo[natom].atomnr = strtol(anr, nullptr, 10);
 727         atoms->pdbinfo[natom].altloc = altloc;
 728         strcpy(atoms->pdbinfo[natom].atomnm, anm_copy);
 729         atoms->pdbinfo[natom].bfac  = strtod(bfac, nullptr);
 730         atoms->pdbinfo[natom].occup = strtod(occup, nullptr);
 731     }
 732     natom++;
 733
 734     return natom;
 735 }
 736
 737 gmx_bool is_hydrogen(const char* nm)
 738 {
 739     char buf[30];
 740
 741     std::strcpy(buf, nm);
 742     trim(buf);
 743
 744     return ((buf[0] == 'H') || ((std::isdigit(buf[0]) != 0) && (buf[1] == 'H')));
 745 }
 746
 747 gmx_bool is_dummymass(const char* nm)
 748 {
 749     char buf[30];
 750
 751     std::strcpy(buf, nm);
 752     trim(buf);
 753
 754     return (buf[0] == 'M') && (std::isdigit(buf[strlen(buf) - 1]) != 0);
 755 }
 756
 757 static void gmx_conect_addline(gmx_conect_t* con, char* line)
 758 {
 759     int n, ai, aj;
 760
 761     std::string form2  = "%*s";
 762     std::string format = form2 + "%d";
 763     if (sscanf(line, format.c_str(), &ai) == 1)
 764     {
 765         do
 766         {
 767             form2 += "%*s";
 768             format = form2 + "%d";
 769             n      = sscanf(line, format.c_str(), &aj);
 770             if (n == 1)
 771             {
 772                 gmx_conect_add(con, ai - 1, aj - 1); /* to prevent duplicated records */
 773             }
 774         } while (n == 1);
 775     }
 776 }
 777
 778 void gmx_conect_dump(FILE* fp, gmx_conect conect)
 779 {
 780     gmx_conect_t* gc = static_cast<gmx_conect_t*>(conect);
 781     int           i;
 782
 783     for (i = 0; (i < gc->nconect); i++)
 784     {
 785         fprintf(fp, "%6s%5d%5d\n", "CONECT", gc->conect[i].ai + 1, gc->conect[i].aj + 1);
 786     }
 787 }
 788
 789 gmx_conect gmx_conect_init()
 790 {
 791     gmx_conect_t* gc;
 792
 793     snew(gc, 1);
 794
 795     return gc;
 796 }
 797
 798 void gmx_conect_done(gmx_conect conect)
 799 {
 800     gmx_conect_t* gc = conect;
 801
 802     sfree(gc->conect);
 803 }
 804
 805 gmx_bool gmx_conect_exist(gmx_conect conect, int ai, int aj)
 806 {
 807     gmx_conect_t* gc = conect;
 808     int           i;
 809
 810     /* if (!gc->bSorted)
 811        sort_conect(gc);*/
 812
 813     for (i = 0; (i < gc->nconect); i++)
 814     {
 815         if (((gc->conect[i].ai == ai) && (gc->conect[i].aj == aj))
 816             || ((gc->conect[i].aj == ai) && (gc->conect[i].ai == aj)))
 817         {
 818             return TRUE;
 819         }
 820     }
 821     return FALSE;
 822 }
 823
 824 void gmx_conect_add(gmx_conect conect, int ai, int aj)
 825 {
 826     gmx_conect_t* gc = static_cast<gmx_conect_t*>(conect);
 827
 828     /* if (!gc->bSorted)
 829        sort_conect(gc);*/
 830
 831     if (!gmx_conect_exist(conect, ai, aj))
 832     {
 833         srenew(gc->conect, ++gc->nconect);
 834         gc->conect[gc->nconect - 1].ai = ai;
 835         gc->conect[gc->nconect - 1].aj = aj;
 836     }
 837 }
 838
 839 int read_pdbfile(FILE*      in,
 840                  char*      title,
 841                  int*       model_nr,
 842                  t_atoms*   atoms,
 843                  t_symtab*  symtab,
 844                  rvec       x[],
 845                  PbcType*   pbcType,
 846                  matrix     box,
 847                  gmx_conect conect)
 848 {
 849     gmx_conect_t* gc = conect;
 850     gmx_bool      bCOMPND;
 851     gmx_bool      bConnWarn = FALSE;
 852     char          line[STRLEN + 1];
 853     char *        c, *d;
 854     int           natom, chainnum;
 855     gmx_bool      bStop = FALSE;
 856
 857     if (pbcType)
 858     {
 859         /* Only assume pbc when there is a CRYST1 entry */
 860         *pbcType = PbcType::No;
 861     }
 862     if (box != nullptr)
 863     {
 864         clear_mat(box);
 865     }
 866
 867     atoms->haveMass    = FALSE;
 868     atoms->haveCharge  = FALSE;
 869     atoms->haveType    = FALSE;
 870     atoms->haveBState  = FALSE;
 871     atoms->havePdbInfo = (atoms->pdbinfo != nullptr);
 872
 873     bCOMPND  = FALSE;
 874     title[0] = '\0';
 875     natom    = 0;
 876     chainnum = 0;
 877     static const gmx::StringToEnumValueConverter<PdbRecordType, enumValueToString, gmx::StringCompareType::Exact, gmx::StripStrings::Yes>
 878             sc_pdbRecordTypeIdentifier;
 879     while (!bStop && (fgets2(line, STRLEN, in) != nullptr))
 880     {
 881         // Convert line to a null-terminated string, then take a substring of at most 6 chars
 882         std::string                  lineAsString(line);
 883         std::optional<PdbRecordType> line_type =
 884                 sc_pdbRecordTypeIdentifier.valueFrom(lineAsString.substr(0, 6));
 885         if (!line_type)
 886         {
 887             // Skip this line because it does not start with a
 888             // recognized leading string describing a PDB record type.
 889             continue;
 890         }
 891
 892         switch (line_type.value())
 893         {
 894             case PdbRecordType::Atom:
 895             case PdbRecordType::Hetatm:
 896                 natom = read_atom(symtab, line, line_type.value(), natom, atoms, x, chainnum);
 897                 break;
 898
 899             case PdbRecordType::Anisou:
 900                 if (atoms->havePdbInfo)
 901                 {
 902                     read_anisou(line, natom, atoms);
 903                 }
 904                 break;
 905
 906             case PdbRecordType::Cryst1: read_cryst1(line, pbcType, box); break;
 907
 908             case PdbRecordType::Title:
 909             case PdbRecordType::Header:
 910                 if (std::strlen(line) > 6)
 911                 {
 912                     c = line + 6;
 913                     /* skip HEADER or TITLE and spaces */
 914                     while (c[0] != ' ')
 915                     {
 916                         c++;
 917                     }
 918                     while (c[0] == ' ')
 919                     {
 920                         c++;
 921                     }
 922                     /* truncate after title */
 923                     d = std::strstr(c, "      ");
 924                     if (d)
 925                     {
 926                         d[0] = '\0';
 927                     }
 928                     if (std::strlen(c) > 0)
 929                     {
 930                         std::strcpy(title, c);
 931                     }
 932                 }
 933                 break;
 934
 935             case PdbRecordType::Compound:
 936                 if ((!std::strstr(line, ": ")) || (std::strstr(line + 6, "MOLECULE:")))
 937                 {
 938                     if (!(c = std::strstr(line + 6, "MOLECULE:")))
 939                     {
 940                         c = line;
 941                     }
 942                     /* skip 'MOLECULE:' and spaces */
 943                     while (c[0] != ' ')
 944                     {
 945                         c++;
 946                     }
 947                     while (c[0] == ' ')
 948                     {
 949                         c++;
 950                     }
 951                     /* truncate after title */
 952                     d = strstr(c, "   ");
 953                     if (d)
 954                     {
 955                         while ((d[-1] == ';') && d > c)
 956                         {
 957                             d--;
 958                         }
 959                         d[0] = '\0';
 960                     }
 961                     if (strlen(c) > 0)
 962                     {
 963                         if (bCOMPND)
 964                         {
 965                             std::strcat(title, "; ");
 966                             std::strcat(title, c);
 967                         }
 968                         else
 969                         {
 970                             std::strcpy(title, c);
 971                         }
 972                     }
 973                     bCOMPND = TRUE;
 974                 }
 975                 break;
 976
 977             case PdbRecordType::Ter: chainnum++; break;
 978
 979             case PdbRecordType::Model:
 980                 if (model_nr)
 981                 {
 982                     sscanf(line, "%*s%d", model_nr);
 983                 }
 984                 break;
 985
 986             case PdbRecordType::EndModel: bStop = TRUE; break;
 987             case PdbRecordType::Conect:
 988                 if (gc)
 989                 {
 990                     gmx_conect_addline(gc, line);
 991                 }
 992                 else if (!bConnWarn)
 993                 {
 994                     fprintf(stderr, "WARNING: all CONECT records are ignored\n");
 995                     bConnWarn = TRUE;
 996                 }
 997                 break;
 998
 999             default: break;
1000         }
1001     }
1002
1003     return natom;
1004 }
1005
1006 void get_pdb_coordnum(FILE* in, int* natoms)
1007 {
1008     char line[STRLEN];
1009
1010     *natoms = 0;
1011     while (fgets2(line, STRLEN, in))
1012     {
1013         if (std::strncmp(line, "ENDMDL", 6) == 0)
1014         {
1015             break;
1016         }
1017         if ((std::strncmp(line, "ATOM  ", 6) == 0) || (std::strncmp(line, "HETATM", 6) == 0))
1018         {
1019             (*natoms)++;
1020         }
1021     }
1022 }
1023
1024 void gmx_pdb_read_conf(const char* infile, t_symtab* symtab, char** name, t_atoms* atoms, rvec x[], PbcType* pbcType, matrix box)
1025 {
1026     FILE* in = gmx_fio_fopen(infile, "r");
1027     char  title[STRLEN];
1028     read_pdbfile(in, title, nullptr, atoms, symtab, x, pbcType, box, nullptr);
1029     if (name != nullptr)
1030     {
1031         *name = gmx_strdup(title);
1032     }
1033     gmx_fio_fclose(in);
1034 }
1035
1036 gmx_conect gmx_conect_generate(const t_topology* top)
1037 {
1038     int        f, i;
1039     gmx_conect gc;
1040
1041     /* Fill the conect records */
1042     gc = gmx_conect_init();
1043
1044     for (f = 0; (f < F_NRE); f++)
1045     {
1046         if (IS_CHEMBOND(f))
1047         {
1048             for (i = 0; (i < top->idef.il[f].nr); i += interaction_function[f].nratoms + 1)
1049             {
1050                 gmx_conect_add(gc, top->idef.il[f].iatoms[i + 1], top->idef.il[f].iatoms[i + 2]);
1051             }
1052         }
1053     }
1054     return gc;
1055 }
1056
1057 int gmx_fprintf_pdb_atomline(FILE*         fp,
1058                              PdbRecordType record,
1059                              int           atom_seq_number,
1060                              const char*   atom_name,
1061                              char          alternate_location,
1062                              const char*   res_name,
1063                              char          chain_id,
1064                              int           res_seq_number,
1065                              char          res_insertion_code,
1066                              real          x,
1067                              real          y,
1068                              real          z,
1069                              real          occupancy,
1070                              real          b_factor,
1071                              const char*   element)
1072 {
1073     char     tmp_atomname[6], tmp_resname[6];
1074     gmx_bool start_name_in_col13;
1075     int      n;
1076
1077     if (record != PdbRecordType::Atom && record != PdbRecordType::Hetatm)
1078     {
1079         gmx_fatal(FARGS, "Can only print PDB atom lines as ATOM or HETATM records");
1080     }
1081
1082     /* Format atom name */
1083     if (atom_name != nullptr)
1084     {
1085         /* If the atom name is an element name with two chars, it should start already in column 13.
1086          * Otherwise it should start in column 14, unless the name length is 4 chars.
1087          */
1088         if ((element != nullptr) && (std::strlen(element) >= 2)
1089             && (gmx_strncasecmp(atom_name, element, 2) == 0))
1090         {
1091             start_name_in_col13 = TRUE;
1092         }
1093         else
1094         {
1095             start_name_in_col13 = (std::strlen(atom_name) >= 4);
1096         }
1097         snprintf(tmp_atomname, sizeof(tmp_atomname), start_name_in_col13 ? "" : " ");
1098         std::strncat(tmp_atomname, atom_name, 4);
1099         tmp_atomname[5] = '\0';
1100     }
1101     else
1102     {
1103         tmp_atomname[0] = '\0';
1104     }
1105
1106     /* Format residue name */
1107     std::strncpy(tmp_resname, (res_name != nullptr) ? res_name : "", 4);
1108     /* Make sure the string is terminated if strlen was > 4 */
1109     tmp_resname[4] = '\0';
1110     /* String is properly terminated, so now we can use strcat. By adding a
1111      * space we can write it right-justified, and if the original name was
1112      * three characters or less there will be a space added on the right side.
1113      */
1114     std::strcat(tmp_resname, " ");
1115
1116     /* Truncate integers so they fit */
1117     atom_seq_number = atom_seq_number % 100000;
1118     res_seq_number  = res_seq_number % 10000;
1119
1120     n = fprintf(fp,
1121                 "%-6s%5d %-4.4s%c%4.4s%c%4d%c   %8.3f%8.3f%8.3f%6.2f%6.2f          %2s\n",
1122                 enumValueToString(record),
1123                 atom_seq_number,
1124                 tmp_atomname,
1125                 alternate_location,
1126                 tmp_resname,
1127                 chain_id,
1128                 res_seq_number,
1129                 res_insertion_code,
1130                 x,
1131                 y,
1132                 z,
1133                 occupancy,
1134                 b_factor,
1135                 (element != nullptr) ? element : "");
1136
1137     return n;
1138 }