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