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