src/gmxlib/pdbio.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * GROningen Mixture of Alchemy and Childrens' Stories
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <ctype.h>
  40
  41 #include "sysstuff.h"
  42 #include "string2.h"
  43 #include "vec.h"
  44 #include "smalloc.h"
  45 #include "typedefs.h"
  46 #include "symtab.h"
  47 #include "pdbio.h"
  48 #include "vec.h"
  49 #include "copyrite.h"
  50 #include "futil.h"
  51 #include "atomprop.h"
  52 #include "physics.h"
  53 #include "pbc.h"
  54 #include "gmxfio.h"
  55
  56 typedef struct {
  57   int ai,aj;
  58 } gmx_conection_t;
  59
  60 typedef struct gmx_conect_t {
  61   int  nconect;
  62   gmx_bool bSorted;
  63   gmx_conection_t *conect;
  64 } gmx_conect_t;
  65
  66 static const char *pdbtp[epdbNR]={
  67   "ATOM  ","HETATM", "ANISOU", "CRYST1",
  68   "COMPND", "MODEL", "ENDMDL", "TER", "HEADER", "TITLE", "REMARK",
  69   "CONECT"
  70 };
  71
  72
  73 /* this is not very good,
  74    but these are only used in gmx_trjconv and gmx_editconv */
  75 static gmx_bool bWideFormat=FALSE;
  76 #define REMARK_SIM_BOX "REMARK    THIS IS A SIMULATION BOX"
  77
  78 void set_pdb_wide_format(gmx_bool bSet)
  79 {
  80   bWideFormat = bSet;
  81 }
  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     temp=name[0];
  91     for(i=1; i<length; i++)
  92       name[i-1]=name[i];
  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                 temp=name[length-1];
 105                 for(i=length-1; i>0; --i)
 106                         name[i]=name[i-1];
 107                 name[0]=temp;
 108         }
 109 }
 110
 111
 112 void gmx_write_pdb_box(FILE *out,int ePBC,matrix box)
 113 {
 114   real alpha,beta,gamma;
 115
 116   if (ePBC == -1)
 117     ePBC = guess_ePBC(box);
 118
 119   if (ePBC == epbcNONE)
 120     return;
 121
 122   if (norm2(box[YY])*norm2(box[ZZ])!=0)
 123     alpha = RAD2DEG*acos(cos_angle_no_table(box[YY],box[ZZ]));
 124   else
 125     alpha = 90;
 126   if (norm2(box[XX])*norm2(box[ZZ])!=0)
 127     beta  = RAD2DEG*acos(cos_angle_no_table(box[XX],box[ZZ]));
 128   else
 129     beta  = 90;
 130   if (norm2(box[XX])*norm2(box[YY])!=0)
 131     gamma = RAD2DEG*acos(cos_angle_no_table(box[XX],box[YY]));
 132   else
 133     gamma = 90;
 134   fprintf(out,"REMARK    THIS IS A SIMULATION BOX\n");
 135   if (ePBC != epbcSCREW) {
 136     fprintf(out,"CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f %-11s%4d\n",
 137             10*norm(box[XX]),10*norm(box[YY]),10*norm(box[ZZ]),
 138             alpha,beta,gamma,"P 1",1);
 139   } else {
 140     /* Double the a-vector length and write the correct space group */
 141     fprintf(out,"CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f %-11s%4d\n",
 142             20*norm(box[XX]),10*norm(box[YY]),10*norm(box[ZZ]),
 143             alpha,beta,gamma,"P 21 1 1",1);
 144
 145   }
 146 }
 147
 148 static void read_cryst1(char *line,int *ePBC,matrix box)
 149 {
 150 #define SG_SIZE 11
 151   char sa[12],sb[12],sc[12],sg[SG_SIZE+1],ident;
 152   double fa,fb,fc,alpha,beta,gamma,cosa,cosb,cosg,sing;
 153   int  syma,symb,symc;
 154   int  ePBC_file;
 155
 156   sscanf(line,"%*s%s%s%s%lf%lf%lf",sa,sb,sc,&alpha,&beta,&gamma);
 157
 158   ePBC_file = -1;
 159   if (strlen(line) >= 55) {
 160     strncpy(sg,line+55,SG_SIZE);
 161     sg[SG_SIZE] = '\0';
 162     ident = ' ';
 163     syma  = 0;
 164     symb  = 0;
 165     symc  = 0;
 166     sscanf(sg,"%c %d %d %d",&ident,&syma,&symb,&symc);
 167     if (ident == 'P' && syma ==  1 && symb <= 1 && symc <= 1) {
 168       fc = strtod(sc,NULL)*0.1;
 169       ePBC_file = (fc > 0 ? epbcXYZ : epbcXY);
 170     }
 171     if (ident == 'P' && syma == 21 && symb == 1 && symc == 1) {
 172       ePBC_file = epbcSCREW;
 173     }
 174   }
 175   if (ePBC) {
 176     *ePBC = ePBC_file;
 177   }
 178
 179   if (box) {
 180     fa = strtod(sa,NULL)*0.1;
 181     fb = strtod(sb,NULL)*0.1;
 182     fc = strtod(sc,NULL)*0.1;
 183     if (ePBC_file == epbcSCREW) {
 184       fa *= 0.5;
 185     }
 186     clear_mat(box);
 187     box[XX][XX] = fa;
 188     if ((alpha!=90.0) || (beta!=90.0) || (gamma!=90.0)) {
 189       if (alpha != 90.0) {
 190         cosa = cos(alpha*DEG2RAD);
 191       } else {
 192         cosa = 0;
 193       }
 194       if (beta != 90.0) {
 195         cosb = cos(beta*DEG2RAD);
 196       } else {
 197         cosb = 0;
 198       }
 199       if (gamma != 90.0) {
 200         cosg = cos(gamma*DEG2RAD);
 201         sing = sin(gamma*DEG2RAD);
 202       } else {
 203         cosg = 0;
 204         sing = 1;
 205       }
 206       box[YY][XX] = fb*cosg;
 207       box[YY][YY] = fb*sing;
 208       box[ZZ][XX] = fc*cosb;
 209       box[ZZ][YY] = fc*(cosa - cosb*cosg)/sing;
 210       box[ZZ][ZZ] = sqrt(fc*fc
 211                          - box[ZZ][XX]*box[ZZ][XX] - box[ZZ][YY]*box[ZZ][YY]);
 212     } else {
 213       box[YY][YY] = fb;
 214       box[ZZ][ZZ] = fc;
 215     }
 216   }
 217 }
 218
 219 void write_pdbfile_indexed(FILE *out,const char *title,
 220                            t_atoms *atoms,rvec x[],
 221                            int ePBC,matrix box,char chainid,
 222                            int model_nr, atom_id nindex, atom_id index[],
 223                            gmx_conect conect, gmx_bool bTerSepChains)
 224 {
 225   gmx_conect_t *gc = (gmx_conect_t *)conect;
 226   char resnm[6],nm[6],pdbform[128],pukestring[100];
 227   atom_id i,ii;
 228   int  resind,resnr,type;
 229   unsigned char resic,ch;
 230   real occup,bfac;
 231   gmx_bool bOccup;
 232   int  nlongname=0;
 233   int  chainnum,lastchainnum;
 234   int  lastresind,lastchainresind;
 235   gmx_residuetype_t rt;
 236   const char *p_restype;
 237   const char *p_lastrestype;
 238
 239   gmx_residuetype_init(&rt);
 240
 241   bromacs(pukestring,99);
 242   fprintf(out,"TITLE     %s\n",(title && title[0])?title:pukestring);
 243   if (bWideFormat) {
 244     fprintf(out,"REMARK    This file does not adhere to the PDB standard\n");
 245     fprintf(out,"REMARK    As a result of, some programs may not like it\n");
 246   }
 247   if (box && ( norm2(box[XX]) || norm2(box[YY]) || norm2(box[ZZ]) ) ) {
 248     gmx_write_pdb_box(out,ePBC,box);
 249   }
 250   if (atoms->pdbinfo) {
 251     /* Check whether any occupancies are set, in that case leave it as is,
 252      * otherwise set them all to one
 253      */
 254     bOccup = TRUE;
 255     for (ii=0; (ii<nindex) && bOccup; ii++) {
 256       i      = index[ii];
 257       bOccup = bOccup && (atoms->pdbinfo[i].occup == 0.0);
 258     }
 259   }
 260   else
 261     bOccup = FALSE;
 262
 263   fprintf(out,"MODEL %8d\n",model_nr>0 ? model_nr : 1);
 264
 265   lastchainresind   = -1;
 266   lastchainnum      = -1;
 267   resind            = -1;
 268   p_restype = NULL;
 269
 270   for (ii=0; ii<nindex; ii++) {
 271     i=index[ii];
 272     lastresind = resind;
 273     resind = atoms->atom[i].resind;
 274     chainnum = atoms->resinfo[resind].chainnum;
 275     p_lastrestype = p_restype;
 276     gmx_residuetype_get_type(rt,*atoms->resinfo[resind].name,&p_restype);
 277
 278     /* Add a TER record if we changed chain, and if either the previous or this chain is protein/DNA/RNA. */
 279     if( bTerSepChains && ii>0 && chainnum != lastchainnum)
 280     {
 281         /* Only add TER if the previous chain contained protein/DNA/RNA. */
 282         if(gmx_residuetype_is_protein(rt,p_lastrestype) || gmx_residuetype_is_dna(rt,p_lastrestype) || gmx_residuetype_is_rna(rt,p_lastrestype))
 283         {
 284             fprintf(out,"TER\n");
 285         }
 286         lastchainnum    = chainnum;
 287         lastchainresind = lastresind;
 288     }
 289
 290     strncpy(resnm,*atoms->resinfo[resind].name,sizeof(resnm)-1);
 291     strncpy(nm,*atoms->atomname[i],sizeof(nm)-1);
 292     /* rename HG12 to 2HG1, etc. */
 293     xlate_atomname_gmx2pdb(nm);
 294     resnr = atoms->resinfo[resind].nr;
 295     resic = atoms->resinfo[resind].ic;
 296     if (chainid!=' ')
 297     {
 298       ch = chainid;
 299     }
 300     else
 301     {
 302       ch = atoms->resinfo[resind].chainid;
 303
 304       if (ch == 0)
 305       {
 306           ch = ' ';
 307       }
 308     }
 309     if (resnr>=10000)
 310       resnr = resnr % 10000;
 311     if (atoms->pdbinfo) {
 312       type  = atoms->pdbinfo[i].type;
 313       occup = bOccup ? 1.0 : atoms->pdbinfo[i].occup;
 314       bfac  = atoms->pdbinfo[i].bfac;
 315     }
 316     else {
 317       type  = 0;
 318       occup = 1.0;
 319       bfac  = 0.0;
 320     }
 321     if (bWideFormat)
 322       strcpy(pdbform,
 323              "%-6s%5u %-4.4s %3.3s %c%4d%c   %10.5f%10.5f%10.5f%8.4f%8.4f    %2s\n");
 324     else {
 325       /* Check whether atomname is an element name */
 326       if ((strlen(nm)<4) && (gmx_strcasecmp(nm,atoms->atom[i].elem) != 0))
 327         strcpy(pdbform,pdbformat);
 328       else {
 329         strcpy(pdbform,pdbformat4);
 330         if (strlen(nm) > 4) {
 331           int maxwln=20;
 332           if (nlongname < maxwln) {
 333             fprintf(stderr,"WARNING: Writing out atom name (%s) longer than 4 characters to .pdb file\n",nm);
 334           } else if (nlongname == maxwln) {
 335             fprintf(stderr,"WARNING: More than %d long atom names, will not write more warnings\n",maxwln);
 336           }
 337           nlongname++;
 338         }
 339       }
 340       strcat(pdbform,"%6.2f%6.2f          %2s\n");
 341     }
 342     fprintf(out,pdbform,pdbtp[type],(i+1)%100000,nm,resnm,ch,resnr,
 343             (resic == '\0') ? ' ' : resic,
 344             10*x[i][XX],10*x[i][YY],10*x[i][ZZ],occup,bfac,atoms->atom[i].elem);
 345     if (atoms->pdbinfo && atoms->pdbinfo[i].bAnisotropic) {
 346       fprintf(out,"ANISOU%5u  %-4.4s%3.3s %c%4d%c %7d%7d%7d%7d%7d%7d\n",
 347               (i+1)%100000,nm,resnm,ch,resnr,
 348               (resic == '\0') ? ' ' : resic,
 349               atoms->pdbinfo[i].uij[0],atoms->pdbinfo[i].uij[1],
 350               atoms->pdbinfo[i].uij[2],atoms->pdbinfo[i].uij[3],
 351               atoms->pdbinfo[i].uij[4],atoms->pdbinfo[i].uij[5]);
 352     }
 353   }
 354
 355   fprintf(out,"TER\n");
 356   fprintf(out,"ENDMDL\n");
 357
 358   if (NULL != gc) {
 359     /* Write conect records */
 360     for(i=0; (i<gc->nconect); i++) {
 361       fprintf(out,"CONECT%5d%5d\n",gc->conect[i].ai+1,gc->conect[i].aj+1);
 362     }
 363   }
 364
 365   gmx_residuetype_destroy(rt);
 366 }
 367
 368 void write_pdbfile(FILE *out,const char *title, t_atoms *atoms,rvec x[],
 369                    int ePBC,matrix box,char chainid,int model_nr,gmx_conect conect,gmx_bool bTerSepChains)
 370 {
 371   atom_id i,*index;
 372
 373   snew(index,atoms->nr);
 374   for(i=0; i<atoms->nr; i++)
 375     index[i]=i;
 376   write_pdbfile_indexed(out,title,atoms,x,ePBC,box,chainid,model_nr,
 377                         atoms->nr,index,conect,bTerSepChains);
 378   sfree(index);
 379 }
 380
 381 int line2type(char *line)
 382 {
 383   int  k;
 384   char type[8];
 385
 386   for(k=0; (k<6); k++)
 387     type[k]=line[k];
 388   type[k]='\0';
 389
 390   for(k=0; (k<epdbNR); k++)
 391     if (strncmp(type,pdbtp[k],strlen(pdbtp[k])) == 0)
 392       break;
 393
 394   return k;
 395 }
 396
 397 static void read_anisou(char line[],int natom,t_atoms *atoms)
 398 {
 399   int  i,j,k,atomnr;
 400   char nc='\0';
 401   char anr[12],anm[12];
 402
 403   /* Skip over type */
 404   j=6;
 405   for(k=0; (k<5); k++,j++) anr[k]=line[j];
 406   anr[k]=nc;
 407   j++;
 408   for(k=0; (k<4); k++,j++) anm[k]=line[j];
 409   anm[k]=nc;
 410   j++;
 411
 412   /* Strip off spaces */
 413   trim(anm);
 414
 415   /* Search backwards for number and name only */
 416   atomnr = strtol(anr, NULL, 10);
 417   for(i=natom-1; (i>=0); i--)
 418     if ((strcmp(anm,*(atoms->atomname[i])) == 0) &&
 419         (atomnr == atoms->pdbinfo[i].atomnr))
 420       break;
 421   if (i < 0)
 422     fprintf(stderr,"Skipping ANISOU record (atom %s %d not found)\n",
 423             anm,atomnr);
 424   else {
 425     if (sscanf(line+29,"%d%d%d%d%d%d",
 426                &atoms->pdbinfo[i].uij[U11],&atoms->pdbinfo[i].uij[U22],
 427                &atoms->pdbinfo[i].uij[U33],&atoms->pdbinfo[i].uij[U12],
 428                &atoms->pdbinfo[i].uij[U13],&atoms->pdbinfo[i].uij[U23])
 429                  == 6) {
 430       atoms->pdbinfo[i].bAnisotropic = TRUE;
 431     }
 432     else {
 433       fprintf(stderr,"Invalid ANISOU record for atom %d\n",i);
 434       atoms->pdbinfo[i].bAnisotropic = FALSE;
 435     }
 436   }
 437 }
 438
 439 void get_pdb_atomnumber(t_atoms *atoms,gmx_atomprop_t aps)
 440 {
 441   int  i,atomnumber,len;
 442   size_t k;
 443   char anm[6],anm_copy[6],*ptr;
 444   char nc='\0';
 445   real eval;
 446
 447   if (!atoms->pdbinfo) {
 448     gmx_incons("Trying to deduce atomnumbers when no pdb information is present");
 449   }
 450   for(i=0; (i<atoms->nr); i++) {
 451     strcpy(anm,atoms->pdbinfo[i].atomnm);
 452     strcpy(anm_copy,atoms->pdbinfo[i].atomnm);
 453     len = strlen(anm);
 454     atomnumber = NOTSET;
 455     if ((anm[0] != ' ') && ((len <=2) || ((len > 2) && !isdigit(anm[2])))) {
 456       anm_copy[2] = nc;
 457       if (gmx_atomprop_query(aps,epropElement,"???",anm_copy,&eval))
 458         atomnumber = gmx_nint(eval);
 459       else {
 460         anm_copy[1] = nc;
 461         if (gmx_atomprop_query(aps,epropElement,"???",anm_copy,&eval))
 462           atomnumber = gmx_nint(eval);
 463       }
 464     }
 465     if (atomnumber == NOTSET) {
 466       k=0;
 467       while ((k < strlen(anm)) && (isspace(anm[k]) || isdigit(anm[k])))
 468         k++;
 469       anm_copy[0] = anm[k];
 470       anm_copy[1] = nc;
 471       if (gmx_atomprop_query(aps,epropElement,"???",anm_copy,&eval))
 472         atomnumber = gmx_nint(eval);
 473     }
 474     atoms->atom[i].atomnumber = atomnumber;
 475     ptr = gmx_atomprop_element(aps,atomnumber);
 476     strncpy(atoms->atom[i].elem,ptr==NULL ? "" : ptr,4);
 477     if (debug)
 478       fprintf(debug,"Atomnumber for atom '%s' is %d\n",anm,atomnumber);
 479   }
 480 }
 481
 482 static int read_atom(t_symtab *symtab,
 483                      char line[],int type,int natom,
 484                      t_atoms *atoms,rvec x[],int chainnum,gmx_bool bChange)
 485 {
 486   t_atom *atomn;
 487   int  j,k;
 488   char nc='\0';
 489   char anr[12],anm[12],anm_copy[12],altloc,resnm[12],rnr[12];
 490   char xc[12],yc[12],zc[12],occup[12],bfac[12];
 491   unsigned char resic;
 492   char chainid;
 493   int  resnr,atomnumber;
 494
 495   if (natom>=atoms->nr)
 496     gmx_fatal(FARGS,"\nFound more atoms (%d) in pdb file than expected (%d)",
 497               natom+1,atoms->nr);
 498
 499   /* Skip over type */
 500   j=6;
 501   for(k=0; (k<5); k++,j++) anr[k]=line[j];
 502   anr[k]=nc;
 503   trim(anr);
 504   j++;
 505   for(k=0; (k<4); k++,j++) anm[k]=line[j];
 506   anm[k]=nc;
 507   strcpy(anm_copy,anm);
 508   rtrim(anm_copy);
 509   atomnumber = NOTSET;
 510   trim(anm);
 511   altloc=line[j];
 512   j++;
 513   for(k=0; (k<4); k++,j++)
 514     resnm[k]=line[j];
 515   resnm[k]=nc;
 516   trim(resnm);
 517
 518   chainid = line[j];
 519   j++;
 520
 521   for(k=0; (k<4); k++,j++) {
 522     rnr[k] = line[j];
 523   }
 524   rnr[k] = nc;
 525   trim(rnr);
 526   resnr = strtol(rnr, NULL, 10);
 527   resic = line[j];
 528   j+=4;
 529
 530   /* X,Y,Z Coordinate */
 531   for(k=0; (k<8); k++,j++) xc[k]=line[j];
 532   xc[k]=nc;
 533   for(k=0; (k<8); k++,j++) yc[k]=line[j];
 534   yc[k]=nc;
 535   for(k=0; (k<8); k++,j++) zc[k]=line[j];
 536   zc[k]=nc;
 537
 538   /* Weight */
 539   for(k=0; (k<6); k++,j++) occup[k]=line[j];
 540   occup[k]=nc;
 541
 542   /* B-Factor */
 543   for(k=0; (k<7); k++,j++) bfac[k]=line[j];
 544   bfac[k]=nc;
 545
 546   if (atoms->atom) {
 547     atomn=&(atoms->atom[natom]);
 548     if ((natom==0) ||
 549         atoms->resinfo[atoms->atom[natom-1].resind].nr != resnr ||
 550         atoms->resinfo[atoms->atom[natom-1].resind].ic != resic ||
 551         (strcmp(*atoms->resinfo[atoms->atom[natom-1].resind].name,resnm) != 0))
 552     {
 553       if (natom == 0) {
 554         atomn->resind = 0;
 555       } else {
 556         atomn->resind = atoms->atom[natom-1].resind + 1;
 557       }
 558       atoms->nres = atomn->resind + 1;
 559       t_atoms_set_resinfo(atoms,natom,symtab,resnm,resnr,resic,chainnum,chainid);
 560     }
 561     else
 562     {
 563       atomn->resind = atoms->atom[natom-1].resind;
 564     }
 565     if (bChange) {
 566       xlate_atomname_pdb2gmx(anm);
 567     }
 568     atoms->atomname[natom]=put_symtab(symtab,anm);
 569     atomn->m = 0.0;
 570     atomn->q = 0.0;
 571     atomn->atomnumber = atomnumber;
 572     atomn->elem[0] = '\0';
 573   }
 574   x[natom][XX]=strtod(xc,NULL)*0.1;
 575   x[natom][YY]=strtod(yc,NULL)*0.1;
 576   x[natom][ZZ]=strtod(zc,NULL)*0.1;
 577   if (atoms->pdbinfo) {
 578     atoms->pdbinfo[natom].type=type;
 579     atoms->pdbinfo[natom].atomnr=strtol(anr, NULL, 10);
 580     atoms->pdbinfo[natom].altloc=altloc;
 581     strcpy(atoms->pdbinfo[natom].atomnm,anm_copy);
 582     atoms->pdbinfo[natom].bfac=strtod(bfac,NULL);
 583     atoms->pdbinfo[natom].occup=strtod(occup,NULL);
 584   }
 585   natom++;
 586
 587   return natom;
 588 }
 589
 590 gmx_bool is_hydrogen(const char *nm)
 591 {
 592   char buf[30];
 593
 594   strcpy(buf,nm);
 595   trim(buf);
 596
 597   if (buf[0] == 'H')
 598     return TRUE;
 599   else if ((isdigit(buf[0])) && (buf[1] == 'H'))
 600     return TRUE;
 601   return FALSE;
 602 }
 603
 604 gmx_bool is_dummymass(const char *nm)
 605 {
 606   char buf[30];
 607
 608   strcpy(buf,nm);
 609   trim(buf);
 610
 611   if ((buf[0] == 'M') && isdigit(buf[strlen(buf)-1]))
 612     return TRUE;
 613
 614   return FALSE;
 615 }
 616
 617 static void gmx_conect_addline(gmx_conect_t *con,char *line)
 618 {
 619   int n,ai,aj;
 620   char format[32],form2[32];
 621
 622   sprintf(form2,"%%*s");
 623   sprintf(format,"%s%%d",form2);
 624   if (sscanf(line,format,&ai) == 1) {
 625     do {
 626       strcat(form2,"%*s");
 627       sprintf(format,"%s%%d",form2);
 628       n = sscanf(line,format,&aj);
 629       if (n == 1) {
 630         srenew(con->conect,++con->nconect);
 631         con->conect[con->nconect-1].ai = ai-1;
 632         con->conect[con->nconect-1].aj = aj-1;
 633       }
 634     } while (n == 1);
 635   }
 636 }
 637
 638 void gmx_conect_dump(FILE *fp,gmx_conect conect)
 639 {
 640   gmx_conect_t *gc = (gmx_conect_t *)conect;
 641   int i;
 642
 643   for(i=0; (i<gc->nconect); i++)
 644     fprintf(fp,"%6s%5d%5d\n","CONECT",
 645             gc->conect[i].ai+1,gc->conect[i].aj+1);
 646 }
 647
 648 gmx_conect gmx_conect_init()
 649 {
 650   gmx_conect_t *gc;
 651
 652   snew(gc,1);
 653
 654   return (gmx_conect) gc;
 655 }
 656
 657 void gmx_conect_done(gmx_conect conect)
 658 {
 659   gmx_conect_t *gc = (gmx_conect_t *)conect;
 660
 661   sfree(gc->conect);
 662 }
 663
 664 gmx_bool gmx_conect_exist(gmx_conect conect,int ai,int aj)
 665 {
 666   gmx_conect_t *gc = (gmx_conect_t *)conect;
 667   int i;
 668
 669   /* if (!gc->bSorted)
 670      sort_conect(gc);*/
 671
 672   for(i=0; (i<gc->nconect); i++)
 673     if (((gc->conect[i].ai == ai) &&
 674          (gc->conect[i].aj == aj)) ||
 675         ((gc->conect[i].aj == ai) &&
 676          (gc->conect[i].ai == aj)))
 677       return TRUE;
 678   return FALSE;
 679 }
 680
 681 void gmx_conect_add(gmx_conect conect,int ai,int aj)
 682 {
 683   gmx_conect_t *gc = (gmx_conect_t *)conect;
 684   int i;
 685
 686   /* if (!gc->bSorted)
 687      sort_conect(gc);*/
 688
 689   if (!gmx_conect_exist(conect,ai,aj)) {
 690     srenew(gc->conect,++gc->nconect);
 691     gc->conect[gc->nconect-1].ai = ai;
 692     gc->conect[gc->nconect-1].aj = aj;
 693   }
 694 }
 695
 696 int read_pdbfile(FILE *in,char *title,int *model_nr,
 697                  t_atoms *atoms,rvec x[],int *ePBC,matrix box,gmx_bool bChange,
 698                  gmx_conect conect)
 699 {
 700     gmx_conect_t *gc = (gmx_conect_t *)conect;
 701     t_symtab symtab;
 702     gmx_bool bCOMPND;
 703     gmx_bool bConnWarn = FALSE;
 704     char line[STRLEN+1];
 705     int  line_type;
 706     char *c,*d;
 707     int  natom,chainnum,nres_ter_prev=0;
 708     char chidmax=' ';
 709     gmx_bool bStop=FALSE;
 710
 711     if (ePBC)
 712     {
 713         /* Only assume pbc when there is a CRYST1 entry */
 714         *ePBC = epbcNONE;
 715     }
 716     if (box != NULL)
 717     {
 718         clear_mat(box);
 719     }
 720
 721     open_symtab(&symtab);
 722
 723     bCOMPND=FALSE;
 724     title[0]='\0';
 725     natom=0;
 726     chainnum=0;
 727     while (!bStop && (fgets2(line,STRLEN,in) != NULL))
 728     {
 729         line_type = line2type(line);
 730
 731         switch(line_type)
 732         {
 733             case epdbATOM:
 734             case epdbHETATM:
 735                 natom = read_atom(&symtab,line,line_type,natom,atoms,x,chainnum,bChange);
 736                 break;
 737
 738             case epdbANISOU:
 739                 if (atoms->pdbinfo)
 740                 {
 741                     read_anisou(line,natom,atoms);
 742                 }
 743                 break;
 744
 745             case epdbCRYST1:
 746                 read_cryst1(line,ePBC,box);
 747                 break;
 748
 749             case epdbTITLE:
 750             case epdbHEADER:
 751                 if (strlen(line) > 6)
 752                 {
 753                     c=line+6;
 754                     /* skip HEADER or TITLE and spaces */
 755                     while (c[0]!=' ') c++;
 756                     while (c[0]==' ') c++;
 757                     /* truncate after title */
 758                     d=strstr(c,"      ");
 759                     if (d)
 760                     {
 761                         d[0]='\0';
 762                     }
 763                     if (strlen(c)>0)
 764                     {
 765                         strcpy(title,c);
 766                     }
 767                 }
 768                 break;
 769
 770             case epdbCOMPND:
 771                 if ((!strstr(line,": ")) || (strstr(line+6,"MOLECULE:")))
 772                 {
 773                     if ( !(c=strstr(line+6,"MOLECULE:")) )
 774                     {
 775                         c=line;
 776                     }
 777                     /* skip 'MOLECULE:' and spaces */
 778                     while (c[0]!=' ') c++;
 779                     while (c[0]==' ') c++;
 780                     /* truncate after title */
 781                     d=strstr(c,"   ");
 782                     if (d)
 783                     {
 784                         while ( (d[-1]==';') && d>c)  d--;
 785                         d[0]='\0';
 786                     }
 787                     if (strlen(c) > 0)
 788                     {
 789                         if (bCOMPND)
 790                         {
 791                             strcat(title,"; ");
 792                             strcat(title,c);
 793                         }
 794                         else
 795                         {
 796                             strcpy(title,c);
 797                         }
 798                     }
 799                     bCOMPND=TRUE;
 800                 }
 801                 break;
 802
 803             case epdbTER:
 804                 chainnum++;
 805                 break;
 806
 807             case epdbMODEL:
 808                 if(model_nr)
 809                 {
 810                     sscanf(line,"%*s%d",model_nr);
 811                 }
 812                 break;
 813
 814             case epdbENDMDL:
 815                 bStop=TRUE;
 816                 break;
 817             case epdbCONECT:
 818                 if (gc)
 819                 {
 820                     gmx_conect_addline(gc,line);
 821                 }
 822                 else if (!bConnWarn)
 823                 {
 824                     fprintf(stderr,"WARNING: all CONECT records are ignored\n");
 825                     bConnWarn = TRUE;
 826                 }
 827                 break;
 828
 829             default:
 830                 break;
 831         }
 832     }
 833
 834     free_symtab(&symtab);
 835     return natom;
 836 }
 837
 838 void get_pdb_coordnum(FILE *in,int *natoms)
 839 {
 840     char line[STRLEN];
 841
 842     *natoms=0;
 843     while (fgets2(line,STRLEN,in))
 844     {
 845         if ( strncmp(line,"ENDMDL",6) == 0 )
 846         {
 847             break;
 848         }
 849         if ((strncmp(line,"ATOM  ",6) == 0) || (strncmp(line,"HETATM",6) == 0))
 850         {
 851             (*natoms)++;
 852         }
 853     }
 854 }
 855
 856 void read_pdb_conf(const char *infile,char *title,
 857                    t_atoms *atoms,rvec x[],int *ePBC,matrix box,gmx_bool bChange,
 858                    gmx_conect conect)
 859 {
 860   FILE *in;
 861
 862   in = gmx_fio_fopen(infile,"r");
 863   read_pdbfile(in,title,NULL,atoms,x,ePBC,box,bChange,conect);
 864   gmx_fio_fclose(in);
 865 }
 866
 867 gmx_conect gmx_conect_generate(t_topology *top)
 868 {
 869   int f,i;
 870   gmx_conect gc;
 871
 872   /* Fill the conect records */
 873   gc  = gmx_conect_init();
 874
 875   for(f=0; (f<F_NRE); f++) {
 876     if (IS_CHEMBOND(f))
 877       for(i=0; (i<top->idef.il[f].nr); i+=interaction_function[f].nratoms+1) {
 878         gmx_conect_add(gc,top->idef.il[f].iatoms[i+1],
 879                        top->idef.il[f].iatoms[i+2]);
 880     }
 881   }
 882   return gc;
 883 }