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