src/gromacs/fileio/pdbio.c

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