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