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