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