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