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