src/tools/editconf.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  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_editconf_c = "$Id$";
  30
  31 #include <math.h>
  32 #include <string.h>
  33 #include <ctype.h>
  34 #include "pdbio.h"
  35 #include "confio.h"
  36 #include "symtab.h"
  37 #include "smalloc.h"
  38 #include "macros.h"
  39 #include "copyrite.h"
  40 #include "statutil.h"
  41 #include "string2.h"
  42 #include "strdb.h"
  43 #include "rdgroup.h"
  44 #include "vec.h"
  45 #include "typedefs.h"
  46 #include "gbutil.h"
  47 #include "strdb.h"
  48 #include "rdgroup.h"
  49 #include "physics.h"
  50 #include "atomprop.h"
  51 #include "addconf.h"
  52 #include "tpxio.h"
  53
  54 typedef struct {
  55   char   sanm[12];
  56   int    natm;
  57   int    nw;
  58   char   anm[6][12];
  59 }  t_simat;
  60
  61 typedef struct {
  62   char     reso[12];
  63   char     resn[12];
  64   int      nsatm;
  65   t_simat sat[3];
  66 } t_simlist;
  67 static char *pdbtp[epdbNR]={"ATOM  ","HETATM"};
  68 static char *pdbformat=
  69 "%6s%5d  %-4.4s%3.3s %c%4d    %8.3f%8.3f%8.3f%6.2f%6.2f\n";
  70
  71 real calc_mass(t_atoms *atoms,bool bGetMass)
  72 {
  73   real tmass;
  74   int i;
  75
  76   tmass = 0;
  77   for(i=0; (i<atoms->nr); i++) {
  78     if (bGetMass)
  79       atoms->atom[i].m = get_mass(*atoms->resname[atoms->atom[i].resnr],
  80                                   *atoms->atomname[i]);
  81     tmass += atoms->atom[i].m;
  82   }
  83
  84   return tmass;
  85 }
  86
  87 void calc_geom(char *indexfn, t_atoms *atoms,
  88                rvec *x, rvec geom_center, rvec min, rvec max)
  89 {
  90   int     isize;
  91   atom_id *index;
  92   char    *grpnames;
  93   int     ii,i,j;
  94
  95   if (indexfn) {
  96     /* Make an index for principal component analysis */
  97     fprintf(stderr,"\nSelect group for selecting system size:\n");
  98     get_index(atoms,indexfn,1,&isize,&index,&grpnames);
  99   }
 100   else {
 101     isize = atoms->nr;
 102     snew(index,isize);
 103     for(i=0; (i<isize); i++) {
 104       index[i]=i;
 105     }
 106   }
 107   clear_rvec(geom_center);
 108   for (j=0; (j<DIM); j++)
 109     min[j]=max[j]=x[0][j];
 110   for (i=0; (i<isize); i++) {
 111     ii = index[i];
 112     rvec_inc(geom_center,x[ii]);
 113     for (j=0; (j<DIM); j++) {
 114       if (x[ii][j] < min[j]) min[j]=x[ii][j];
 115       if (x[ii][j] > max[j]) max[j]=x[ii][j];
 116     }
 117   }
 118   svmul(1.0/isize,geom_center,geom_center);
 119   sfree(index);
 120 }
 121
 122 void center_conf(int natom, rvec *x, rvec center, rvec geom_cent)
 123 {
 124   int       i;
 125   rvec shift;
 126
 127   rvec_sub(center,geom_cent,shift);
 128
 129   printf("shift     : %6.3f %6.3f %6.3f\n",
 130          shift[XX],shift[YY],shift[ZZ]);
 131
 132   for (i=0; (i<natom); i++)
 133     rvec_inc(x[i], shift);
 134 }
 135
 136 void scale_conf(int natom,rvec x[],matrix box,rvec scale)
 137 {
 138   int i,j;
 139
 140   for(i=0; (i<natom); i++) {
 141     for (j=0; (j<DIM); j++)
 142       x[i][j] *= scale[j];
 143   }
 144   for (j=0; (j<DIM); j++)
 145     box[j][j] *= scale[j];
 146 }
 147
 148 void rm_gropbc(t_atoms *atoms,rvec x[],matrix box)
 149 {
 150   real dist;
 151   int  n,d;
 152
 153   /* check periodic boundary */
 154   for(d=0;(d<DIM);d++) {
 155     for(n=1;(n<atoms->nr);n++) {
 156       dist = x[n][d]-x[n-1][d];
 157       if ( fabs(dist) > 0.9 * box[d][d]  ) {
 158         if ( dist >  0 )
 159           x[n][d]-=box[d][d];
 160         else
 161           x[n][d]+=box[d][d];
 162       }
 163     }
 164   }
 165 }
 166
 167 void read_bfac(char *fn, int *n_bfac, double **bfac_val, int **bfac_nr)
 168 {
 169   int  i;
 170   char **bfac_lines;
 171
 172   *n_bfac = get_lines(fn, &bfac_lines);
 173   snew(*bfac_val, *n_bfac);
 174   snew(*bfac_nr, *n_bfac);
 175   fprintf(stderr, "Reading %d B-factors from %s\n",*n_bfac,fn);
 176   for(i=0; (i<*n_bfac); i++) {
 177     /*fprintf(stderr, "Line %d: %s",i,bfac_lines[i]);*/
 178     sscanf(bfac_lines[i],"%d %lf",&(*bfac_nr)[i],&(*bfac_val)[i]);
 179     /*fprintf(stderr," nr %d val %g\n",(*bfac_nr)[i],(*bfac_val)[i]);*/
 180   }
 181
 182 }
 183
 184 void set_pdb_conf_bfac(int natoms,int nres,t_atoms *atoms,rvec x[],
 185                        matrix box,int n_bfac,double *bfac,int *bfac_nr,
 186                        bool peratom, bool bLegend)
 187 {
 188   FILE *out;
 189   real bfac_min,bfac_max;
 190   int  i,n;
 191   bool found;
 192   char buf[120];
 193
 194   bfac_max=-1e10;
 195   bfac_min=1e10;
 196   for(i=0; (i<n_bfac); i++) {
 197     if (bfac_nr[i]-1>=atoms->nres)
 198       peratom=TRUE;
 199     if ((bfac_nr[i]-1<0) || (bfac_nr[i]-1>=atoms->nr))
 200       fatal_error(0,"Index of B-Factor %d is out of range: %d (%g)",
 201                   i+1,bfac_nr[i],bfac[i]);
 202     if (bfac[i] > bfac_max)
 203       bfac_max = bfac[i];
 204     if (bfac[i] < bfac_min)
 205       bfac_min = bfac[i];
 206   }
 207   while ( (bfac_max > 99.99) || (bfac_min < -99.99) ) {
 208     fprintf(stderr,"Range of values for B-factors too large (min %g, max %g) "
 209             "will scale down a factor 10\n",bfac_min,bfac_max);
 210     for(i=0; (i<n_bfac); i++)
 211       bfac[i] /= 10;
 212     bfac_max /= 10;
 213     bfac_min /= 10;
 214   }
 215   while ( (abs(bfac_max) < 0.5) && (abs(bfac_min) < 0.5) ) {
 216     fprintf(stderr,"Range of values for B-factors too small (min %g, max %g) "
 217             "will scale up a factor 10\n",bfac_min,bfac_max);
 218     for(i=0; (i<n_bfac); i++)
 219       bfac[i] *= 10;
 220     bfac_max *= 10;
 221     bfac_min *= 10;
 222   }
 223
 224   for(i=0; (i<natoms); i++)
 225     atoms->pdbinfo[i].bfac=0;
 226
 227   if (!peratom) {
 228     fprintf(stderr,"Will attach %d B-factors to %d residues\n",
 229             n_bfac,nres);
 230     for(i=0; (i<n_bfac); i++) {
 231       found=FALSE;
 232       for(n=0; (n<natoms); n++)
 233         if ( bfac_nr[i] == (atoms->atom[n].resnr+1) ) {
 234           atoms->pdbinfo[n].bfac=bfac[i];
 235           found=TRUE;
 236         }
 237       if (!found) {
 238         sprintf(buf,"Residue nr %d not found\n",bfac_nr[i]);
 239         warning(buf);
 240       }
 241     }
 242   } else {
 243     fprintf(stderr,"Will attach %d B-factors to %d atoms\n",n_bfac,natoms);
 244     for(i=0; (i<n_bfac); i++) {
 245       atoms->pdbinfo[bfac_nr[i]-1].bfac=bfac[i];
 246     }
 247   }
 248 }
 249
 250 void pdb_legend(FILE *out,int natoms,int nres,t_atoms *atoms,rvec x[])
 251 {
 252   real bfac_min,bfac_max,xmin,ymin,zmin;
 253   int  i;
 254   char buf[256];
 255
 256   bfac_max=-1e10;
 257   bfac_min=1e10;
 258   xmin = 1e10;
 259   ymin = 1e10;
 260   zmin = 1e10;
 261   for (i=0; (i<natoms); i++) {
 262     xmin     = min(xmin,x[i][XX]);
 263     ymin     = min(ymin,x[i][YY]);
 264     zmin     = min(zmin,x[i][ZZ]);
 265     bfac_min = min(bfac_min,atoms->pdbinfo[i].bfac);
 266     bfac_max = max(bfac_max,atoms->pdbinfo[i].bfac);
 267   }
 268   fprintf(stderr,"B-factors range from %g to %g\n",bfac_min,bfac_max);
 269   sprintf(buf,"%s","LEG");
 270   buf[3]='\0';
 271   for (i=1; (i<12); i++) {
 272     fprintf(out,pdbformat,
 273             "ATOM  ",natoms+1+i,"CA",buf,' ',nres+1,
 274             (xmin+(i*0.12))*10,ymin*10,zmin*10,1.0,
 275             bfac_min+ ((i-1.0)*(bfac_max-bfac_min)/10) );
 276   }
 277 }
 278
 279 int main(int argc, char *argv[])
 280 {
 281   static char *desc[] = {
 282     "editconf converts generic structure format to [TT].gro[tt] or",
 283     "[TT].pdb[tt].[PAR]",
 284     "A number of options is present to modify the coordinates",
 285     "and box. [TT]-d[tt], [TT]-dc[tt] and [TT]-box[tt] modify the box and",
 286     "center the coordinates relative to the new box.",
 287     "[TT]-dc[tt] takes precedent over [TT]-d[tt]. [TT]-box[tt]",
 288     "takes precedent over [TT]-dc[tt] and [TT]-d[tt].[PAR]",
 289     "[TT]-rotate[tt] rotates the coordinates and velocities.",
 290     "[TT]-princ[tt] aligns the principal axes of the system along the",
 291     "coordinate axes, this may allow you to decrease the box volume,",
 292     "but beware that molecules can rotate significantly in a nanosecond.[PAR]",
 293     "Scaling is applied before any of the other operations are",
 294     "performed. Boxes can be scaled to give a certain density (option",
 295     "[TT]-density[tt]). A special feature of the scaling option, when the",
 296     "factor -1 is given in one dimension, one obtains a mirror image,",
 297     "mirrored in one of the plains, when one uses -1 in three dimensions",
 298     "a point-mirror image is obtained.[PAR]",
 299     "Groups are selected after all operations have been applied.[PAR]",
 300     "Periodicity can be removed in a crude manner.",
 301     "It is important that the box sizes at the bottom of your input file",
 302     "are correct when the periodicity is to be removed.[PAR]",
 303     "The program can optionally rotate the solute molecule to align the",
 304     "molecule along its principal axes ([TT]-rotate[tt])[PAR]",
 305     "When writing [TT].pdb[tt] files, B-factors can be",
 306     "added with the [TT]-bf[tt] option. B-factors are read",
 307     "from a file with with following format: first line states number of",
 308     "entries in the file, next lines state an index",
 309     "followed by a B-factor. The B-factors will be attached per residue",
 310     "unless an index is larger than the number of residues or unless the",
 311     "[TT]-atom[tt] option is set. Obviously, any type of numeric data can",
 312     "be added instead of B-factors. [TT]-legend[tt] will produce",
 313     "a row of CA atoms with B-factors ranging from the minimum to the",
 314     "maximum value found, effectively making a legend for viewing.[PAR]",
 315     "Finally with option [TT]-label[tt] editconf can add a chain identifier",
 316     "to a pdb file, which can be useful for analysis with e.g. rasmol."
 317   };
 318   static char *bugs[] = {
 319     "For complex molecules, the periodicity removal routine may break down, "
 320     "in that case you can use trjconv"
 321   };
 322   static real dist   = 0.0,rbox=0.0;
 323   static bool bNDEF=FALSE,bRMPBC=FALSE,bCenter=FALSE;
 324   static bool peratom=FALSE,bLegend=FALSE,bOrient=FALSE;
 325   static rvec scale={1.0,1.0,1.0},newbox={0.0,0.0,0.0};
 326   static real rho=1000.0;
 327   static rvec center={0.0,0.0,0.0},rotangles={0.0,0.0,0.0};
 328   static char *label="A";
 329   t_pargs pa[] = {
 330     { "-ndef",   FALSE, etBOOL, {&bNDEF},
 331       "Choose output from default index groups" },
 332     { "-d",      FALSE, etREAL, {&dist},
 333       "Distance between the solute and the rectangular box" },
 334     { "-dc",     FALSE, etREAL, {&dist},
 335       "Distance between the solute and the cubic box" },
 336     { "-box",    FALSE, etRVEC, {&newbox}, "Size of box" },
 337     { "-c",      FALSE, etBOOL, {&bCenter},
 338       "Center molecule in box (implied by -d -dc -box)" },
 339     { "-center", FALSE, etRVEC, {&center}, "Coordinates of geometrical center"},
 340     { "-rotate", FALSE, etRVEC, {rotangles},
 341       "Rotation around the X, Y and Z axes in degrees" },
 342     { "-princ",  FALSE, etBOOL, {&bOrient}, "Orient molecule(s) along their principal axes" },
 343     { "-scale",  FALSE, etRVEC, {&scale}, "Scaling factor" },
 344     { "-density",FALSE, etREAL, {&rho},
 345       "Density (g/l) of the output box achieved by scaling" },
 346     { "-pbc",    FALSE, etBOOL, {&bRMPBC},
 347       "Remove the periodicity (make molecule whole again)" },
 348     { "-atom",   FALSE, etBOOL, {&peratom}, "Force B-factor attachment per atom" },
 349     { "-legend", FALSE, etBOOL, {&bLegend}, "Make B-factor legend" },
 350     { "-label",  FALSE, etSTR,  {&label},   "Add chain label for all residues" }
 351   };
 352 #define NPA asize(pa)
 353
 354   FILE      *out;
 355   char      *infile,*outfile,title[STRLEN];
 356   int       outftp,natom,i,j,n_bfac;
 357   double    *bfac=NULL;
 358   int       *bfac_nr=NULL;
 359   t_atoms   atoms;
 360   char      *groupnames;
 361   int       isize;
 362   atom_id   *index;
 363   rvec      *x,*v,gc,min,max,size;
 364   matrix    box;
 365   bool      bSetSize,bCubic,bDist,bSetCenter;
 366   bool      bHaveV,bScale,bRho,bRotate,bCalcGeom;
 367   real      xs,ys,zs,xcent,ycent,zcent,d;
 368   t_filenm fnm[] = {
 369     { efSTX, "-f", NULL, ffREAD },
 370     { efNDX, "-n", NULL, ffOPTRD },
 371     { efSTO, NULL, NULL, ffWRITE },
 372     { efDAT, "-bf", "bfact", ffOPTRD }
 373   };
 374 #define NFILE asize(fnm)
 375
 376   CopyRight(stderr,argv[0]);
 377   parse_common_args(&argc,argv,0,FALSE,NFILE,fnm,NPA,pa,
 378                     asize(desc),desc,asize(bugs),bugs);
 379
 380   bSetSize  = opt2parg_bSet("-box" ,NPA,pa);
 381   bSetCenter= opt2parg_bSet("-center" ,NPA,pa);
 382   bCubic    = opt2parg_bSet("-dc",NPA,pa);
 383   bDist     = opt2parg_bSet("-d" ,NPA,pa) || bCubic;
 384   bCenter   = bCenter || bDist || bSetCenter || bSetSize;
 385   bScale    = opt2parg_bSet("-scale" ,NPA,pa);
 386   bRho      = opt2parg_bSet("-density",NPA,pa);
 387   bRotate   = opt2parg_bSet("-rotate",NPA,pa);
 388   if (bScale && bRho)
 389     fprintf(stderr,"WARNING: setting -density overrides -scale");
 390   bScale    = bScale || bRho;
 391   bCalcGeom = bCenter || bRotate || bOrient || bScale;
 392
 393   infile  = ftp2fn(efSTX,NFILE,fnm);
 394   outfile = ftp2fn(efSTO,NFILE,fnm);
 395   outftp  = fn2ftp(outfile);
 396
 397   get_stx_coordnum(infile,&natom);
 398   init_t_atoms(&atoms,natom,TRUE);
 399   snew(x,natom);
 400   snew(v,natom);
 401   read_stx_conf(infile,title,&atoms,x,v,box);
 402   printf("Read %d atoms\n",atoms.nr);
 403
 404   bHaveV=FALSE;
 405   for (i=0; (i<natom) && !bHaveV; i++)
 406     for (j=0; (j<DIM) && !bHaveV; j++)
 407       bHaveV=bHaveV || (v[i][j]!=0);
 408   printf("%selocities found\n",bHaveV?"V":"No v");
 409
 410   /* remove pbc */
 411   if (bRMPBC)
 412     rm_gropbc(&atoms,x,box);
 413
 414   if (bCalcGeom) {
 415     calc_geom(ftp2fn_null(efNDX,NFILE,fnm),&atoms,x, gc, min, max);
 416     rvec_sub(max, min, size);
 417     printf("size      : %6.3f %6.3f %6.3f\n", size[XX], size[YY], size[ZZ]);
 418     printf("center    : %6.3f %6.3f %6.3f\n", gc[XX], gc[YY], gc[ZZ]);
 419     printf("box       : %6.3f %6.3f %6.3f  (%.3f nm^3)\n",
 420            box[XX][XX], box[YY][YY], box[ZZ][ZZ], det(box));
 421   }
 422
 423   if (bOrient)
 424     /* Orient the principal axes along the coordinate axes */
 425     orient_mol(&atoms,ftp2fn_null(efNDX,NFILE,fnm),x,bHaveV ? v : NULL);
 426
 427   if ( bScale ) {
 428     /* scale coordinates and box */
 429     if (bRho) {
 430       /* Compute scaling constant */
 431       real vol,mass,dens;
 432
 433       vol = det(box);
 434       mass = calc_mass(&atoms,!fn2bTPX(infile));
 435       dens = (mass*AMU)/(vol*NANO*NANO*NANO);
 436       fprintf(stderr,"Volume  of input %g (nm3)\n",vol);
 437       fprintf(stderr,"Mass    of input %g (a.m.u.)\n",mass);
 438       fprintf(stderr,"Density of input %g (g/l)\n",dens);
 439       if (vol==0.0)
 440         fatal_error(0,"Cannot scale density with zero box\n");
 441       if (mass==0.0)
 442         fatal_error(0,"Cannot scale density with zero mass\n");
 443
 444       scale[XX] = scale[YY] = scale[ZZ] = pow(dens/rho,1.0/3.0);
 445       fprintf(stderr,"Scaling all box edges by %g\n",scale[XX]);
 446     }
 447     scale_conf(atoms.nr,x,box,scale);
 448   }
 449
 450   if (bRotate) {
 451     /* Rotate */
 452     printf("Rotating %g, %g, %g degrees around the X, Y and Z axis respectively\n",rotangles[XX],rotangles[YY],rotangles[ZZ]);
 453     for(i=0; i<DIM; i++)
 454       rotangles[i] *= DEG2RAD;
 455     rotate_conf(natom,x,v,rotangles[XX],rotangles[YY],rotangles[ZZ]);
 456   }
 457
 458   if (bCalcGeom) {
 459     /* recalc geometrical center and max and min coordinates and size */
 460     calc_geom(ftp2fn_null(efNDX,NFILE,fnm),&atoms,x, gc, min, max);
 461     rvec_sub(max, min, size);
 462     if (bScale || bOrient || bRotate)
 463       printf("new size  : %6.3f %6.3f %6.3f\n",size[XX],size[YY],size[ZZ]);
 464   }
 465
 466   /* calculate new boxsize */
 467   if (bDist)
 468     for (i=0; (i<DIM); i++)
 469       box[i][i]=size[i]+2*dist;
 470   if (bSetSize)
 471     for (i=0; (i<DIM); i++)
 472       box[i][i]=newbox[i];
 473   if (bCubic) {
 474     d=box[XX][XX];
 475     for (i=1; (i<DIM); i++)
 476       if (box[i][i]>d) d=box[i][i];
 477     for (i=0; (i<DIM); i++)
 478       box[i][i]=d;
 479   }
 480   /* calculate new coords for geometrical center */
 481   if (!bSetCenter)
 482     for (i=0; (i<DIM); i++)
 483       center[i]=box[i][i]/2;
 484
 485   /* center molecule on 'center' */
 486   if (bCenter)
 487     center_conf(natom,x,center,gc);
 488
 489   /* print some */
 490   if (bCenter || bScale || bOrient || bRotate) {
 491     calc_geom(NULL,&atoms,x, gc, min, max);
 492     printf("new center: %6.3f %6.3f %6.3f\n",
 493            gc[XX],gc[YY],gc[ZZ]);
 494   }
 495   if ( bOrient || bScale || bDist || bSetSize )
 496     printf("new box   : %6.3f %6.3f %6.3f  (%.3f nm^3)\n",
 497            box[XX][XX], box[YY][YY], box[ZZ][ZZ], det(box));
 498
 499   if (opt2bSet("-n",NFILE,fnm) || bNDEF) {
 500     fprintf(stderr,"\nSelect a group for output:\n");
 501     get_index(&atoms,opt2fn_null("-n",NFILE,fnm),
 502               1,&isize,&index,&groupnames);
 503     if (opt2bSet("-bf",NFILE,fnm))
 504       fatal_error(0,"combination not implemented: -bf -n  or -bf -ndef");
 505     else
 506       write_sto_conf_indexed(outfile,title,&atoms,x,bHaveV?v:NULL,box,
 507                              isize,index);
 508   }
 509   else {
 510     if (outftp != efPDB) {
 511       write_sto_conf(outfile,title,&atoms,x,bHaveV?v:NULL,box);
 512     } else {
 513       out=ffopen(outfile,"w");
 514       if (opt2bSet("-bf",NFILE,fnm)) {
 515         read_bfac(opt2fn("-bf",NFILE,fnm),&n_bfac,&bfac,&bfac_nr);
 516         set_pdb_conf_bfac(atoms.nr,atoms.nres,&atoms,x,box,
 517                           n_bfac,bfac,bfac_nr,peratom,bLegend);
 518       }
 519       if (opt2parg_bSet("-label",NPA,pa)) {
 520         for(i=0; (i<atoms.nr); i++)
 521           atoms.atom[i].chain=label[0];
 522       }
 523       write_pdbfile(out,title,&atoms,x,box,0,!bLegend);
 524       if (bLegend)
 525         pdb_legend(out,atoms.nr,atoms.nres,&atoms,x);
 526       fclose(out);
 527     }
 528   }
 529
 530   thanx(stdout);
 531
 532   return 0;
 533 }
 534