src/gromacs/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,get_pdbformat());
 328       else {
 329         strcpy(pdbform,get_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
 366 void write_pdbfile(FILE *out,const char *title, t_atoms *atoms,rvec x[],
 367                    int ePBC,matrix box,char chainid,int model_nr,gmx_conect conect,gmx_bool bTerSepChains)
 368 {
 369   atom_id i,*index;
 370
 371   snew(index,atoms->nr);
 372   for(i=0; i<atoms->nr; i++)
 373     index[i]=i;
 374   write_pdbfile_indexed(out,title,atoms,x,ePBC,box,chainid,model_nr,
 375                         atoms->nr,index,conect,bTerSepChains);
 376   sfree(index);
 377 }
 378
 379 int line2type(char *line)
 380 {
 381   int  k;
 382   char type[8];
 383
 384   for(k=0; (k<6); k++)
 385     type[k]=line[k];
 386   type[k]='\0';
 387
 388   for(k=0; (k<epdbNR); k++)
 389     if (strncmp(type,pdbtp[k],strlen(pdbtp[k])) == 0)
 390       break;
 391
 392   return k;
 393 }
 394
 395 static void read_anisou(char line[],int natom,t_atoms *atoms)
 396 {
 397   int  i,j,k,atomnr;
 398   char nc='\0';
 399   char anr[12],anm[12];
 400
 401   /* Skip over type */
 402   j=6;
 403   for(k=0; (k<5); k++,j++) anr[k]=line[j];
 404   anr[k]=nc;
 405   j++;
 406   for(k=0; (k<4); k++,j++) anm[k]=line[j];
 407   anm[k]=nc;
 408   j++;
 409
 410   /* Strip off spaces */
 411   trim(anm);
 412
 413   /* Search backwards for number and name only */
 414   atomnr = strtol(anr, NULL, 10);
 415   for(i=natom-1; (i>=0); i--)
 416     if ((strcmp(anm,*(atoms->atomname[i])) == 0) &&
 417         (atomnr == atoms->pdbinfo[i].atomnr))
 418       break;
 419   if (i < 0)
 420     fprintf(stderr,"Skipping ANISOU record (atom %s %d not found)\n",
 421             anm,atomnr);
 422   else {
 423     if (sscanf(line+29,"%d%d%d%d%d%d",
 424                &atoms->pdbinfo[i].uij[U11],&atoms->pdbinfo[i].uij[U22],
 425                &atoms->pdbinfo[i].uij[U33],&atoms->pdbinfo[i].uij[U12],
 426                &atoms->pdbinfo[i].uij[U13],&atoms->pdbinfo[i].uij[U23])
 427                  == 6) {
 428       atoms->pdbinfo[i].bAnisotropic = TRUE;
 429     }
 430     else {
 431       fprintf(stderr,"Invalid ANISOU record for atom %d\n",i);
 432       atoms->pdbinfo[i].bAnisotropic = FALSE;
 433     }
 434   }
 435 }
 436
 437 void get_pdb_atomnumber(t_atoms *atoms,gmx_atomprop_t aps)
 438 {
 439   int  i,atomnumber,len;
 440   size_t k;
 441   char anm[6],anm_copy[6],*ptr;
 442   char nc='\0';
 443   real eval;
 444
 445   if (!atoms->pdbinfo) {
 446     gmx_incons("Trying to deduce atomnumbers when no pdb information is present");
 447   }
 448   for(i=0; (i<atoms->nr); i++) {
 449     strcpy(anm,atoms->pdbinfo[i].atomnm);
 450     strcpy(anm_copy,atoms->pdbinfo[i].atomnm);
 451     len = strlen(anm);
 452     atomnumber = NOTSET;
 453     if ((anm[0] != ' ') && ((len <=2) || ((len > 2) && !isdigit(anm[2])))) {
 454       anm_copy[2] = nc;
 455       if (gmx_atomprop_query(aps,epropElement,"???",anm_copy,&eval))
 456         atomnumber = gmx_nint(eval);
 457       else {
 458         anm_copy[1] = nc;
 459         if (gmx_atomprop_query(aps,epropElement,"???",anm_copy,&eval))
 460           atomnumber = gmx_nint(eval);
 461       }
 462     }
 463     if (atomnumber == NOTSET) {
 464       k=0;
 465       while ((k < strlen(anm)) && (isspace(anm[k]) || isdigit(anm[k])))
 466         k++;
 467       anm_copy[0] = anm[k];
 468       anm_copy[1] = nc;
 469       if (gmx_atomprop_query(aps,epropElement,"???",anm_copy,&eval))
 470         atomnumber = gmx_nint(eval);
 471     }
 472     atoms->atom[i].atomnumber = atomnumber;
 473     ptr = gmx_atomprop_element(aps,atomnumber);
 474     strncpy(atoms->atom[i].elem,ptr==NULL ? "" : ptr,4);
 475     if (debug)
 476       fprintf(debug,"Atomnumber for atom '%s' is %d\n",anm,atomnumber);
 477   }
 478 }
 479
 480 static int read_atom(t_symtab *symtab,
 481                      char line[],int type,int natom,
 482                      t_atoms *atoms,rvec x[],int chainnum,gmx_bool bChange)
 483 {
 484   t_atom *atomn;
 485   int  j,k;
 486   char nc='\0';
 487   char anr[12],anm[12],anm_copy[12],altloc,resnm[12],rnr[12];
 488   char xc[12],yc[12],zc[12],occup[12],bfac[12];
 489   unsigned char resic;
 490   char chainid;
 491   int  resnr,atomnumber;
 492
 493   if (natom>=atoms->nr)
 494     gmx_fatal(FARGS,"\nFound more atoms (%d) in pdb file than expected (%d)",
 495               natom+1,atoms->nr);
 496
 497   /* Skip over type */
 498   j=6;
 499   for(k=0; (k<5); k++,j++) anr[k]=line[j];
 500   anr[k]=nc;
 501   trim(anr);
 502   j++;
 503   for(k=0; (k<4); k++,j++) anm[k]=line[j];
 504   anm[k]=nc;
 505   strcpy(anm_copy,anm);
 506   rtrim(anm_copy);
 507   atomnumber = NOTSET;
 508   trim(anm);
 509   altloc=line[j];
 510   j++;
 511   for(k=0; (k<4); k++,j++)
 512     resnm[k]=line[j];
 513   resnm[k]=nc;
 514   trim(resnm);
 515
 516   chainid = line[j];
 517   j++;
 518
 519   for(k=0; (k<4); k++,j++) {
 520     rnr[k] = line[j];
 521   }
 522   rnr[k] = nc;
 523   trim(rnr);
 524   resnr = strtol(rnr, NULL, 10);
 525   resic = line[j];
 526   j+=4;
 527
 528   /* X,Y,Z Coordinate */
 529   for(k=0; (k<8); k++,j++) xc[k]=line[j];
 530   xc[k]=nc;
 531   for(k=0; (k<8); k++,j++) yc[k]=line[j];
 532   yc[k]=nc;
 533   for(k=0; (k<8); k++,j++) zc[k]=line[j];
 534   zc[k]=nc;
 535
 536   /* Weight */
 537   for(k=0; (k<6); k++,j++) occup[k]=line[j];
 538   occup[k]=nc;
 539
 540   /* B-Factor */
 541   for(k=0; (k<7); k++,j++) bfac[k]=line[j];
 542   bfac[k]=nc;
 543
 544   if (atoms->atom) {
 545     atomn=&(atoms->atom[natom]);
 546     if ((natom==0) ||
 547         atoms->resinfo[atoms->atom[natom-1].resind].nr != resnr ||
 548         atoms->resinfo[atoms->atom[natom-1].resind].ic != resic ||
 549         (strcmp(*atoms->resinfo[atoms->atom[natom-1].resind].name,resnm) != 0))
 550     {
 551       if (natom == 0) {
 552         atomn->resind = 0;
 553       } else {
 554         atomn->resind = atoms->atom[natom-1].resind + 1;
 555       }
 556       atoms->nres = atomn->resind + 1;
 557       t_atoms_set_resinfo(atoms,natom,symtab,resnm,resnr,resic,chainnum,chainid);
 558     }
 559     else
 560     {
 561       atomn->resind = atoms->atom[natom-1].resind;
 562     }
 563     if (bChange) {
 564       xlate_atomname_pdb2gmx(anm);
 565     }
 566     atoms->atomname[natom]=put_symtab(symtab,anm);
 567     atomn->m = 0.0;
 568     atomn->q = 0.0;
 569     atomn->atomnumber = atomnumber;
 570     atomn->elem[0] = '\0';
 571   }
 572   x[natom][XX]=strtod(xc,NULL)*0.1;
 573   x[natom][YY]=strtod(yc,NULL)*0.1;
 574   x[natom][ZZ]=strtod(zc,NULL)*0.1;
 575   if (atoms->pdbinfo) {
 576     atoms->pdbinfo[natom].type=type;
 577     atoms->pdbinfo[natom].atomnr=strtol(anr, NULL, 10);
 578     atoms->pdbinfo[natom].altloc=altloc;
 579     strcpy(atoms->pdbinfo[natom].atomnm,anm_copy);
 580     atoms->pdbinfo[natom].bfac=strtod(bfac,NULL);
 581     atoms->pdbinfo[natom].occup=strtod(occup,NULL);
 582   }
 583   natom++;
 584
 585   return natom;
 586 }
 587
 588 gmx_bool is_hydrogen(const char *nm)
 589 {
 590   char buf[30];
 591
 592   strcpy(buf,nm);
 593   trim(buf);
 594
 595   if (buf[0] == 'H')
 596     return TRUE;
 597   else if ((isdigit(buf[0])) && (buf[1] == 'H'))
 598     return TRUE;
 599   return FALSE;
 600 }
 601
 602 gmx_bool is_dummymass(const char *nm)
 603 {
 604   char buf[30];
 605
 606   strcpy(buf,nm);
 607   trim(buf);
 608
 609   if ((buf[0] == 'M') && isdigit(buf[strlen(buf)-1]))
 610     return TRUE;
 611
 612   return FALSE;
 613 }
 614
 615 static void gmx_conect_addline(gmx_conect_t *con,char *line)
 616 {
 617   int n,ai,aj;
 618   char format[32],form2[32];
 619
 620   sprintf(form2,"%%*s");
 621   sprintf(format,"%s%%d",form2);
 622   if (sscanf(line,format,&ai) == 1) {
 623     do {
 624       strcat(form2,"%*s");
 625       sprintf(format,"%s%%d",form2);
 626       n = sscanf(line,format,&aj);
 627       if (n == 1) {
 628         srenew(con->conect,++con->nconect);
 629         con->conect[con->nconect-1].ai = ai-1;
 630         con->conect[con->nconect-1].aj = aj-1;
 631       }
 632     } while (n == 1);
 633   }
 634 }
 635
 636 void gmx_conect_dump(FILE *fp,gmx_conect conect)
 637 {
 638   gmx_conect_t *gc = (gmx_conect_t *)conect;
 639   int i;
 640
 641   for(i=0; (i<gc->nconect); i++)
 642     fprintf(fp,"%6s%5d%5d\n","CONECT",
 643             gc->conect[i].ai+1,gc->conect[i].aj+1);
 644 }
 645
 646 gmx_conect gmx_conect_init()
 647 {
 648   gmx_conect_t *gc;
 649
 650   snew(gc,1);
 651
 652   return (gmx_conect) gc;
 653 }
 654
 655 void gmx_conect_done(gmx_conect conect)
 656 {
 657   gmx_conect_t *gc = (gmx_conect_t *)conect;
 658
 659   sfree(gc->conect);
 660 }
 661
 662 gmx_bool gmx_conect_exist(gmx_conect conect,int ai,int aj)
 663 {
 664   gmx_conect_t *gc = (gmx_conect_t *)conect;
 665   int i;
 666
 667   /* if (!gc->bSorted)
 668      sort_conect(gc);*/
 669
 670   for(i=0; (i<gc->nconect); i++)
 671     if (((gc->conect[i].ai == ai) &&
 672          (gc->conect[i].aj == aj)) ||
 673         ((gc->conect[i].aj == ai) &&
 674          (gc->conect[i].ai == aj)))
 675       return TRUE;
 676   return FALSE;
 677 }
 678
 679 void gmx_conect_add(gmx_conect conect,int ai,int aj)
 680 {
 681   gmx_conect_t *gc = (gmx_conect_t *)conect;
 682   int i;
 683
 684   /* if (!gc->bSorted)
 685      sort_conect(gc);*/
 686
 687   if (!gmx_conect_exist(conect,ai,aj)) {
 688     srenew(gc->conect,++gc->nconect);
 689     gc->conect[gc->nconect-1].ai = ai;
 690     gc->conect[gc->nconect-1].aj = aj;
 691   }
 692 }
 693
 694 int read_pdbfile(FILE *in,char *title,int *model_nr,
 695                  t_atoms *atoms,rvec x[],int *ePBC,matrix box,gmx_bool bChange,
 696                  gmx_conect conect)
 697 {
 698     gmx_conect_t *gc = (gmx_conect_t *)conect;
 699     t_symtab symtab;
 700     gmx_bool bCOMPND;
 701     gmx_bool bConnWarn = FALSE;
 702     char line[STRLEN+1];
 703     int  line_type;
 704     char *c,*d;
 705     int  natom,chainnum,nres_ter_prev=0;
 706     char chidmax=' ';
 707     gmx_bool bStop=FALSE;
 708
 709     if (ePBC)
 710     {
 711         /* Only assume pbc when there is a CRYST1 entry */
 712         *ePBC = epbcNONE;
 713     }
 714     if (box != NULL)
 715     {
 716         clear_mat(box);
 717     }
 718
 719     open_symtab(&symtab);
 720
 721     bCOMPND=FALSE;
 722     title[0]='\0';
 723     natom=0;
 724     chainnum=0;
 725     while (!bStop && (fgets2(line,STRLEN,in) != NULL))
 726     {
 727         line_type = line2type(line);
 728
 729         switch(line_type)
 730         {
 731             case epdbATOM:
 732             case epdbHETATM:
 733                 natom = read_atom(&symtab,line,line_type,natom,atoms,x,chainnum,bChange);
 734                 break;
 735
 736             case epdbANISOU:
 737                 if (atoms->pdbinfo)
 738                 {
 739                     read_anisou(line,natom,atoms);
 740                 }
 741                 break;
 742
 743             case epdbCRYST1:
 744                 read_cryst1(line,ePBC,box);
 745                 break;
 746
 747             case epdbTITLE:
 748             case epdbHEADER:
 749                 if (strlen(line) > 6)
 750                 {
 751                     c=line+6;
 752                     /* skip HEADER or TITLE and spaces */
 753                     while (c[0]!=' ') c++;
 754                     while (c[0]==' ') c++;
 755                     /* truncate after title */
 756                     d=strstr(c,"      ");
 757                     if (d)
 758                     {
 759                         d[0]='\0';
 760                     }
 761                     if (strlen(c)>0)
 762                     {
 763                         strcpy(title,c);
 764                     }
 765                 }
 766                 break;
 767
 768             case epdbCOMPND:
 769                 if ((!strstr(line,": ")) || (strstr(line+6,"MOLECULE:")))
 770                 {
 771                     if ( !(c=strstr(line+6,"MOLECULE:")) )
 772                     {
 773                         c=line;
 774                     }
 775                     /* skip 'MOLECULE:' and spaces */
 776                     while (c[0]!=' ') c++;
 777                     while (c[0]==' ') c++;
 778                     /* truncate after title */
 779                     d=strstr(c,"   ");
 780                     if (d)
 781                     {
 782                         while ( (d[-1]==';') && d>c)  d--;
 783                         d[0]='\0';
 784                     }
 785                     if (strlen(c) > 0)
 786                     {
 787                         if (bCOMPND)
 788                         {
 789                             strcat(title,"; ");
 790                             strcat(title,c);
 791                         }
 792                         else
 793                         {
 794                             strcpy(title,c);
 795                         }
 796                     }
 797                     bCOMPND=TRUE;
 798                 }
 799                 break;
 800
 801             case epdbTER:
 802                 chainnum++;
 803                 break;
 804
 805             case epdbMODEL:
 806                 if(model_nr)
 807                 {
 808                     sscanf(line,"%*s%d",model_nr);
 809                 }
 810                 break;
 811
 812             case epdbENDMDL:
 813                 bStop=TRUE;
 814                 break;
 815             case epdbCONECT:
 816                 if (gc)
 817                 {
 818                     gmx_conect_addline(gc,line);
 819                 }
 820                 else if (!bConnWarn)
 821                 {
 822                     fprintf(stderr,"WARNING: all CONECT records are ignored\n");
 823                     bConnWarn = TRUE;
 824                 }
 825                 break;
 826
 827             default:
 828                 break;
 829         }
 830     }
 831
 832     free_symtab(&symtab);
 833     return natom;
 834 }
 835
 836 void get_pdb_coordnum(FILE *in,int *natoms)
 837 {
 838     char line[STRLEN];
 839
 840     *natoms=0;
 841     while (fgets2(line,STRLEN,in))
 842     {
 843         if ( strncmp(line,"ENDMDL",6) == 0 )
 844         {
 845             break;
 846         }
 847         if ((strncmp(line,"ATOM  ",6) == 0) || (strncmp(line,"HETATM",6) == 0))
 848         {
 849             (*natoms)++;
 850         }
 851     }
 852 }
 853
 854 void read_pdb_conf(const char *infile,char *title,
 855                    t_atoms *atoms,rvec x[],int *ePBC,matrix box,gmx_bool bChange,
 856                    gmx_conect conect)
 857 {
 858   FILE *in;
 859
 860   in = gmx_fio_fopen(infile,"r");
 861   read_pdbfile(in,title,NULL,atoms,x,ePBC,box,bChange,conect);
 862   gmx_fio_fclose(in);
 863 }
 864
 865 gmx_conect gmx_conect_generate(t_topology *top)
 866 {
 867   int f,i;
 868   gmx_conect gc;
 869
 870   /* Fill the conect records */
 871   gc  = gmx_conect_init();
 872
 873   for(f=0; (f<F_NRE); f++) {
 874     if (IS_CHEMBOND(f))
 875       for(i=0; (i<top->idef.il[f].nr); i+=interaction_function[f].nratoms+1) {
 876         gmx_conect_add(gc,top->idef.il[f].iatoms[i+1],
 877                        top->idef.il[f].iatoms[i+2]);
 878     }
 879   }
 880   return gc;
 881 }
 882
 883 const char* get_pdbformat()
 884 {
 885     static const char *pdbformat ="%-6s%5u  %-4.4s%3.3s %c%4d%c   %8.3f%8.3f%8.3f";
 886     return pdbformat;
 887 }
 888
 889 const char* get_pdbformat4()
 890 {
 891     static const char *pdbformat4="%-6s%5u %-4.4s %3.3s %c%4d%c   %8.3f%8.3f%8.3f";
 892     return pdbformat4;
 893 }