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