src/gromacs/fileio/pdbio.cpp

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