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