src/gromacs/fileio/pdbio.cpp

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