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