src/gmxlib/pdbio.c

   1 /*
   2  * $Id$
   3  *
   4  *       This source code is part of
   5  *
   6  *        G   R   O   M   A   C   S
   7  *
   8  * GROningen MAchine for Chemical Simulations
   9  *
  10  *               VERSION 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * Green Red Orange Magenta Azure Cyan Skyblue
  28  */
  29 static char *SRCID_pdbio_c = "$Id$";
  30
  31 #include "sysstuff.h"
  32 #include "string2.h"
  33 #include "vec.h"
  34 #include "smalloc.h"
  35 #include "typedefs.h"
  36 #include "symtab.h"
  37 #include "assert.h"
  38 #include "pdbio.h"
  39 #include "vec.h"
  40 #include "copyrite.h"
  41 #include "futil.h"
  42
  43 static char *pdbtp[epdbNR]={
  44   "ATOM  ","HETATM", "ANISOU", "CRYST1",
  45   "COMPND", "ENDMDL", "TER", "HEADER", "TITLE", "REMARK"
  46 };
  47
  48 static char *pdbformat ="%6s%5u  %-4.4s%3.3s %c%4d    %8.3f%8.3f%8.3f%6.2f%6.2f\n";
  49 static char *pdbformat4="%6s%5u %-4.4s %3.3s %c%4d    %8.3f%8.3f%8.3f%6.2f%6.2f\n";
  50 static bool bTER=FALSE;
  51 #define REMARK_SIM_BOX "REMARK    THIS IS A SIMULATION BOX"
  52
  53
  54 void pdb_use_ter(bool bSet)
  55 {
  56   bTER=bSet;
  57 }
  58
  59 static void gromacs_name(char *name)
  60 {
  61   int i,length;
  62   char temp;
  63
  64   length=strlen(name);
  65   if (isdigit(name[0])) {
  66     temp=name[0];
  67     for(i=1; i<length; i++)
  68       name[i-1]=name[i];
  69     name[length-1]=temp;
  70   }
  71   if(strcmp(name,"OXT")==0)
  72     strcpy(name,"O2");
  73 }
  74
  75 void write_pdbfile_indexed(FILE *out,char *title,
  76                            t_atoms *atoms,rvec x[],matrix box,char chain,
  77                            bool bEndmodel,
  78                            atom_id nindex, atom_id index[])
  79 {
  80   char resnm[6],nm[6],ch,*pdbform;
  81   atom_id i,ii;
  82   int  resnr,type;
  83   real occup,bfac;
  84
  85   fprintf(out,"HEADER    %s\n",(title && title[0])?title:bromacs());
  86   if (box) {
  87     fprintf(out,"REMARK    THIS IS A SIMULATION BOX\n");
  88     fprintf(out,"CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f P 1           1\n",
  89             10*box[XX][XX],10*box[YY][YY],10*box[ZZ][ZZ],90.0,90.0,90.0);
  90   }
  91   for (ii=0; ii<nindex; ii++) {
  92     i=index[ii];
  93     resnr=atoms->atom[i].resnr;
  94     strcpy(resnm,*atoms->resname[resnr]);
  95     strcpy(nm,*atoms->atomname[i]);
  96     resnr++;
  97     if (resnr>=10000)
  98       resnr = resnr % 10000;
  99     if (chain)
 100       ch=chain;
 101     else
 102       if (atoms->atom[i].chain)
 103         ch=atoms->atom[i].chain;
 104       else
 105           ch=' ';
 106     if (atoms->pdbinfo) {
 107       type=atoms->pdbinfo[i].type;
 108       occup=atoms->pdbinfo[i].occup;
 109       bfac=atoms->pdbinfo[i].bfac;
 110       if (occup==0.0)
 111         occup=1.0;
 112     }
 113     else {
 114       type=0;
 115       occup=1.0;
 116       bfac=0.0;
 117     }
 118     if (strlen(nm)==4)
 119       pdbform=pdbformat4;
 120     else
 121       pdbform=pdbformat;
 122     fprintf(out,pdbform,pdbtp[type],(i+1)%100000,nm,resnm,ch,resnr,
 123     10*x[i][XX],10*x[i][YY],10*x[i][ZZ],occup,bfac);
 124   }
 125
 126   fprintf(out,"TER\n");
 127   if (bEndmodel && !bTER)
 128     fprintf(out,"ENDMDL\n");
 129 }
 130
 131 void write_pdbfile(FILE *out,char *title,
 132                    t_atoms *atoms,rvec x[],matrix box,char chain,
 133                    bool bEndmodel)
 134 {
 135   atom_id i,*index;
 136
 137   snew(index,atoms->nr);
 138   for(i=0; i<atoms->nr; i++)
 139     index[i]=i;
 140   write_pdbfile_indexed(out,title,atoms,x,box,chain,bEndmodel,
 141                         atoms->nr,index);
 142   sfree(index);
 143 }
 144
 145 void hwrite_pdb_conf_indexed(FILE *out,char *title,
 146                              t_atoms *atoms,rvec x[],matrix box,
 147                              atom_id nindex,atom_id index[])
 148 {
 149   write_pdbfile_indexed(out,title,atoms,x,box,0,TRUE,nindex,index);
 150 }
 151
 152 void write_pdb_confs(char *outfile,t_atoms **atoms,rvec *x[],int number)
 153 {
 154   FILE *out;
 155   int n;
 156   char chain,str[STRLEN];
 157
 158   out=ffopen(outfile,"w");
 159
 160   for(n=0;(n<number);n++) {
 161     chain='A'+n;
 162     fprintf(stderr,"writing chain %c\n",chain);
 163     sprintf(str,"Chain %c",chain);
 164     write_pdbfile(out,str,atoms[n],x[n],NULL,chain,(n==number-1));
 165   }
 166
 167   ffclose(out);
 168 }
 169
 170 int line2type(char *line)
 171 {
 172   int  k;
 173   char type[8];
 174
 175   for(k=0; (k<6); k++)
 176     type[k]=line[k];
 177   type[k]='\0';
 178
 179   for(k=0; (k<epdbNR); k++)
 180     if (strncmp(type,pdbtp[k],strlen(pdbtp[k])) == 0)
 181       break;
 182
 183   return k;
 184 }
 185
 186 static void read_anisou(char line[],int natom,t_atoms *atoms)
 187 {
 188   int  i,j,k,atomnr;
 189   char nc='\0';
 190   char anr[12],anm[12];
 191
 192   /* Skip over type */
 193   j=6;
 194   for(k=0; (k<5); k++,j++) anr[k]=line[j];
 195   anr[k]=nc;
 196   j++;
 197   for(k=0; (k<4); k++,j++) anm[k]=line[j];
 198   anm[k]=nc;
 199   j++;
 200
 201   /* Strip off spaces */
 202   trim(anm);
 203
 204   /* Search backwards for number and name only */
 205   atomnr = atoi(anr);
 206   for(i=natom-1; (i>=0); i--)
 207     if ((strcmp(anm,*(atoms->atomname[i])) == 0) &&
 208         (atomnr == atoms->pdbinfo[i].atomnr))
 209       break;
 210   if (i < 0)
 211     fprintf(stderr,"Skipping ANISOU record (atom %s %d not found)\n",
 212             anm,atomnr);
 213   else {
 214     if (sscanf(line+29,"%d%d%d%d%d%d",
 215                &atoms->pdbinfo[i].uij[U11],&atoms->pdbinfo[i].uij[U22],
 216                &atoms->pdbinfo[i].uij[U33],&atoms->pdbinfo[i].uij[U12],
 217                &atoms->pdbinfo[i].uij[U13],&atoms->pdbinfo[i].uij[U23])
 218                  == 6) {
 219       atoms->pdbinfo[i].bAnisotropic = TRUE;
 220     }
 221     else {
 222       fprintf(stderr,"Invalid ANISOU record for atom %d\n",i);
 223       atoms->pdbinfo[i].bAnisotropic = FALSE;
 224     }
 225   }
 226 }
 227
 228 static int read_atom(t_symtab *symtab,char line[],int type,int natom,
 229                      t_atoms *atoms,rvec x[],bool bChange)
 230 {
 231   t_atom *atomn;
 232   int  j,k;
 233   char nc='\0';
 234   char anr[12],anm[12],altloc,resnm[12],chain[12],resnr[12];
 235   char xc[12],yc[12],zc[12],occup[12],bfac[12],pdbresnr[12];
 236   static char oldresnm[12],oldresnr[12];
 237   int  newres;
 238
 239   if (natom>=atoms->nr)
 240     fatal_error(0,"\nFound more atoms (%d) in pdb file than expected (%d)",
 241                 natom+1,atoms->nr);
 242
 243   /* Skip over type */
 244   j=6;
 245   for(k=0; (k<5); k++,j++) anr[k]=line[j];
 246   anr[k]=nc;
 247   trim(anr);
 248   j++;
 249   for(k=0; (k<4); k++,j++) anm[k]=line[j];
 250   anm[k]=nc;
 251   trim(anm);
 252   altloc=line[j];
 253   j++;
 254   for(k=0; (k<4); k++,j++)
 255     resnm[k]=line[j];
 256   resnm[k]=nc;
 257   trim(resnm);
 258
 259   for(k=0; (k<1); k++,j++)
 260     chain[k]=line[j];
 261   chain[k]=nc;
 262
 263   for(k=0; (k<4); k++,j++) {
 264     resnr[k]=line[j];
 265     pdbresnr[k]=line[j];
 266   }
 267   resnr[k]=nc;
 268   trim(resnr);
 269   pdbresnr[k]=line[j];
 270   pdbresnr[k+1]=nc;
 271   trim(pdbresnr);
 272   j+=4;
 273
 274   /* X,Y,Z Coordinate */
 275   for(k=0; (k<8); k++,j++) xc[k]=line[j];
 276   xc[k]=nc;
 277   for(k=0; (k<8); k++,j++) yc[k]=line[j];
 278   yc[k]=nc;
 279   for(k=0; (k<8); k++,j++) zc[k]=line[j];
 280   zc[k]=nc;
 281
 282   /* Weight */
 283   for(k=0; (k<6); k++,j++) occup[k]=line[j];
 284   occup[k]=nc;
 285
 286   /* B-Factor */
 287   for(k=0; (k<6); k++,j++) bfac[k]=line[j];
 288   bfac[k]=nc;
 289
 290   atomn=&(atoms->atom[natom]);
 291   if ((natom==0) || (strcmp(oldresnr,pdbresnr)!=0) ||
 292       (strcmp(oldresnm,resnm)!=0)) {
 293     strcpy(oldresnr,pdbresnr);
 294     strcpy(oldresnm,resnm);
 295     if (natom==0)
 296       newres=0;
 297     else
 298       newres=atoms->atom[natom-1].resnr+1;
 299     atoms->nres=newres+1;
 300     atoms->resname[newres]=put_symtab(symtab,resnm);
 301   }
 302   else
 303     newres=atoms->atom[natom-1].resnr;
 304   if (bChange)
 305     gromacs_name(anm);
 306   atoms->atomname[natom]=put_symtab(symtab,anm);
 307   atomn->chain=chain[0];
 308   atomn->resnr=newres;
 309   x[natom][XX]=atof(xc)*0.1;
 310   x[natom][YY]=atof(yc)*0.1;
 311   x[natom][ZZ]=atof(zc)*0.1;
 312   if (atoms->pdbinfo) {
 313     atoms->pdbinfo[natom].type=type;
 314     atoms->pdbinfo[natom].atomnr=atoi(anr);
 315     atoms->pdbinfo[natom].altloc=altloc;
 316     strcpy(atoms->pdbinfo[natom].pdbresnr,pdbresnr);
 317     atoms->pdbinfo[natom].bfac=atof(bfac);
 318     atoms->pdbinfo[natom].occup=atof(occup);
 319   }
 320   atomn->m = 0.0;
 321   atomn->q = 0.0;
 322   natom++;
 323
 324   return natom;
 325 }
 326
 327 bool is_hydrogen(char *nm)
 328 {
 329   char buf[30];
 330
 331   strcpy(buf,nm);
 332   trim(buf);
 333
 334   if (buf[0] == 'H')
 335     return TRUE;
 336   else if ((isdigit(buf[0])) && (buf[1] == 'H'))
 337     return TRUE;
 338   return FALSE;
 339 }
 340
 341 bool is_dummymass(char *nm)
 342 {
 343   char buf[30];
 344
 345   strcpy(buf,nm);
 346   trim(buf);
 347
 348   if ((buf[0] == 'M') && isdigit(buf[strlen(buf)-1]))
 349     return TRUE;
 350
 351   return FALSE;
 352 }
 353
 354 int read_pdbfile(FILE *in,char *title,
 355                  t_atoms *atoms,rvec x[],matrix box,bool bChange)
 356 {
 357   static t_symtab symtab;
 358   static bool bFirst=TRUE;
 359   bool bCOMPND;
 360   char line[STRLEN+1];
 361   char xc[12],yc[12],zc[12];
 362   double xa,ya,za;
 363   int  line_type;
 364   char *c,*d;
 365   int  natom;
 366   bool bStop=FALSE;
 367
 368   if (box != NULL)
 369     clear_mat(box);
 370
 371   if (bFirst) {
 372     open_symtab(&symtab);
 373     bFirst=FALSE;
 374   }
 375
 376   bCOMPND=FALSE;
 377   title[0]='\0';
 378   natom=0;
 379   while (!bStop && (fgets2(line,STRLEN,in) != NULL)) {
 380     line_type = line2type(line);
 381
 382     switch(line_type) {
 383     case epdbATOM:
 384     case epdbHETATM:
 385       natom = read_atom(&symtab,line,line_type,natom,atoms,x,bChange);
 386       break;
 387
 388     case epdbANISOU:
 389       if (atoms->pdbinfo != NULL)
 390         read_anisou(line,natom,atoms);
 391       break;
 392
 393     case epdbCRYST1:
 394       if (box) {
 395         sscanf(line,"%*s%s%s%s%lf%lf%lf",xc,yc,zc,&xa,&ya,&za);
 396         if (xa==90 || ya==90 || za==90) {
 397           box[XX][XX] = atof(xc)*0.1;
 398           box[YY][YY] = atof(yc)*0.1;
 399           box[ZZ][ZZ] = atof(zc)*0.1;
 400         } else
 401           clear_mat(box);
 402       }
 403       break;
 404
 405     case epdbTITLE:
 406     case epdbHEADER:
 407       c=line+6;
 408       /* skip HEADER or TITLE and spaces */
 409       while (c && (c[0]!=' ')) c++;
 410       while (c && (c[0]==' ')) c++;
 411       /* truncate after title */
 412       d=strstr(c,"      ");
 413       if (d) {
 414         d[0]='\0';
 415       }
 416       if (strlen(c)>0)
 417         strcpy(title,c);
 418       break;
 419
 420     case epdbCOMPND:
 421       if ((!strstr(line,": ")) || (strstr(line+6,"MOLECULE:"))) {
 422         if ( !(c=strstr(line+6,"MOLECULE:")) )
 423           c=line;
 424         /* skip 'MOLECULE:' and spaces */
 425         while (c && (c[0]!=' ')) c++;
 426         while (c && (c[0]==' ')) c++;
 427         /* truncate after title */
 428         d=strstr(c,"   ");
 429         if (d) {
 430           while ( (d[-1]==';') && d>c)  d--;
 431           d[0]='\0';
 432         }
 433         if (strlen(c) > 0) {
 434           if (bCOMPND) {
 435             strcat(title,"; ");
 436             strcat(title,c);
 437           } else
 438             strcpy(title,c);
 439         }
 440         bCOMPND=TRUE;
 441       }
 442       break;
 443
 444     case epdbTER:
 445       if (bTER)
 446         bStop=TRUE;
 447       break;
 448     case epdbENDMDL:
 449       bStop=TRUE;
 450       break;
 451     default:
 452       break;
 453     }
 454   }
 455
 456   return natom;
 457 }
 458
 459 void get_pdb_coordnum(FILE *in,int *natoms)
 460 {
 461   char line[STRLEN];
 462
 463   *natoms=0;
 464   while (fgets2(line,STRLEN,in)) {
 465     if ( ( bTER && (strncmp(line,"TER",3) == 0)) ||
 466          (!bTER && (strncmp(line,"ENDMDL",6) == 0)) )
 467       break;
 468     if ((strncmp(line,"ATOM  ",6) == 0) || (strncmp(line,"HETATM",6) == 0))
 469       (*natoms)++;
 470   }
 471 }
 472
 473 void read_pdb_conf(char *infile,char *title,
 474                    t_atoms *atoms,rvec x[],matrix box,bool bChange)
 475 {
 476   FILE      *in;
 477
 478   in = ffopen(infile,"r");
 479   read_pdbfile(in,title,atoms,x,box,bChange);
 480   ffclose(in);
 481 }