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