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