src/tools/gmx_sdf.c

   1 /*
   2  *
   3  *                This source code is NOT part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  * This program is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU General Public License
  11  * as published by the Free Software Foundation; either version 2
  12  * of the License, or (at your option) any later version.
  13  *
  14  * And Hey:
  15  * Gyas ROwers Mature At Cryogenic Speed
  16  */
  17
  18
  19 #include <math.h>
  20 #include <ctype.h>
  21 #include "string2.h"
  22 #include "sysstuff.h"
  23 #include "typedefs.h"
  24 #include "macros.h"
  25 #include "vec.h"
  26 #include "pbc.h"
  27 #include "rmpbc.h"
  28 #include "copyrite.h"
  29 #include "futil.h"
  30 #include "statutil.h"
  31 #include "tpxio.h"
  32 #include "rdgroup.h"
  33 #include "smalloc.h"
  34 #include "nrnb.h"
  35 #include "gstat.h"
  36 #include "matio.h"
  37 #include "gmx_fatal.h"
  38 #include "gmx_ana.h"
  39
  40
  41 #define G_REF1      0
  42 #define G_REF2      1
  43 #define G_REF3      2
  44 #define G_SDF       3
  45 #define NDX_REF1    4
  46 #define NDX_REF2    5
  47 #define NDX_REF3    6
  48 #define G_REFMOL    7
  49
  50 static void i_write(FILE *output, int value)
  51 {
  52   if(fwrite(&value,sizeof(int),1,output) != 1)
  53   {
  54     gmx_fatal(FARGS,"Error writing to output file");
  55   }
  56 }
  57
  58
  59 static void f_write(FILE *output,float value)
  60 {
  61   if(fwrite(&value,sizeof(float),1,output) != 1)
  62   {
  63     gmx_fatal(FARGS,"Error writing to output file");
  64   }
  65 }
  66
  67
  68 static void do_sdf(const char *fnNDX,const char *fnTPS,const char *fnTRX,
  69                    const char *fnSDF, const char *fnREF, bool bRef,
  70                    rvec cutoff, real binwidth, int mode, rvec triangle,
  71                    rvec dtri, const output_env_t oenv)
  72 {
  73   FILE       *fp;
  74   int        status;
  75   int        ng,natoms,i,j,k,l,X,Y,Z,lc,dest;
  76   ivec       nbin;
  77   int        ***count;
  78   /* real       ***sdf; */
  79   real       sdf,min_sdf=1e10,max_sdf=0;
  80   char       **grpname;
  81   int        *isize;
  82   int        isize_cg=0;
  83   int        isize_ref=3;
  84   int        ref_resind[3]={0};
  85   int        nrefmol=0,refc=0;
  86   atom_id    **index;
  87   atom_id    *index_cg=NULL;
  88   atom_id    *index_ref=NULL;
  89   real       t,boxmin,hbox,normfac;
  90   real       invbinw;
  91   rvec       tri_upper,tri_lower;
  92   rvec       *x,xcog,dx,*x_i1,xi,*x_refmol;
  93   matrix     box;
  94   matrix     rot; /* rotation matrix := unit vectors for the molecule frame */
  95   rvec       k_mol,i1_mol,i2_mol,dx_mol;
  96   real       delta;
  97   atom_id    ix,jx;
  98   t_topology top;
  99   int        ePBC=-1;
 100   t_pbc      pbc;
 101   bool       bTop=FALSE,bRefDone=FALSE,bInGroup=FALSE;
 102   char       title[STRLEN];
 103
 104
 105   /* Read Topology */
 106   if (fnTPS) {
 107     bTop=read_tps_conf(fnTPS,title,&top,&ePBC,&x,NULL,box,TRUE);
 108   }
 109
 110
 111
 112   if ( !bTop ) {
 113     fprintf(stderr,"\nNeed tpr-file to make a reference structure.\n");
 114     fprintf(stderr,"Option -r will be ignored!\n");
 115     bRef = FALSE;
 116   }
 117
 118
 119   /* Allocate memory for 4 groups, 3 dummy groups and a group for the ref
 120 structure if needed */
 121   ng = 4;
 122   snew(grpname,ng);
 123   /* the dummy groups are used to dynamically store triples of atoms */
 124   /* for molecular coordinate systems */
 125   if ( bRef )
 126     {
 127       snew(isize,ng+4);
 128       snew(index,ng+4);
 129     }
 130   else
 131     {
 132       snew(isize,ng+3);
 133       snew(index,ng+3);
 134     }
 135
 136
 137   /* Read the index groups */
 138   fprintf(stderr,"\nSelect the 3 reference groups and the SDF group:\n");
 139   if (fnTPS)
 140     get_index(&top.atoms,fnNDX,ng,isize,index,grpname);
 141   else
 142     rd_index(fnNDX,ng,isize,index,grpname);
 143
 144
 145   isize[NDX_REF1]=isize[G_REF1];
 146   for (i=NDX_REF1; i<=NDX_REF3; i++)
 147     snew(index[i],isize[NDX_REF1]);
 148
 149
 150   /* Read first frame and check it */
 151   natoms=read_first_x(oenv,&status,fnTRX,&t,&x,box);
 152   if ( !natoms )
 153     gmx_fatal(FARGS,"Could not read coordinates from statusfile!\n");
 154
 155
 156   /* check with topology */
 157   if (fnTPS)
 158     if ( natoms > top.atoms.nr )
 159       gmx_fatal(FARGS,
 160                 "Trajectory (%d atoms) does not match topology (%d atoms)!\n",
 161                 natoms,top.atoms.nr);
 162
 163
 164   /* check with index groups */
 165   for (i=0; i<ng; i++)
 166     for (j=0; j<isize[i]; j++)
 167       if ( index[i][j] >= natoms )
 168         gmx_fatal(FARGS,"Atom index (%d) in index group %s (%d atoms) larger "
 169                     "than number of atoms in trajectory (%d atoms)!\n",
 170                     index[i][j],grpname[i],isize[i],natoms);
 171
 172
 173   /* check reference groups */
 174   if ( mode == 1 )
 175     {
 176       if ( isize[G_REF1] != isize[G_REF2] || isize[G_REF1] != isize[G_REF3] ||
 177            isize[G_REF2] != isize[G_REF3] )
 178         gmx_fatal(FARGS,"For single particle SDF, all reference groups"
 179                     "must have the same size.\n");
 180
 181
 182       /* for single particle SDF dynamic triples are not needed */
 183       /* so we build them right here */
 184
 185
 186       /* copy all triples from G_REFx to NDX_REFx */
 187       for (i=0; i<isize[G_REF1]; i++)
 188         {
 189           /* check if all three atoms come from the same molecule */
 190           for (j=G_REF1; j<=G_REF3; j++)
 191             ref_resind[j] = top.atoms.atom[index[j][i]].resind;
 192
 193
 194           if ( ref_resind[G_REF1] != ref_resind[G_REF2] ||
 195                  ref_resind[G_REF2] != ref_resind[G_REF3] ||
 196                  ref_resind[G_REF3] != ref_resind[G_REF1] )
 197               {
 198                 fprintf(stderr,"\nWarning: reference triple (%d) will be skipped.\n",i);
 199                 fprintf(stderr,  "         resnr[1]: %d, resnr[2]: %d, resnr[3]: %d\n",
 200                         ref_resind[G_REF1],ref_resind[G_REF2], ref_resind[G_REF3]);
 201                 isize[NDX_REF1]--;
 202                 for (j=NDX_REF1; j<=NDX_REF3; j++)
 203                   srenew(index[j],isize[NDX_REF1]);
 204                 continue;
 205               }
 206           else
 207             /* check if all entries are unique*/
 208             if ( index[G_REF1][i] == index[G_REF2][i] ||
 209                  index[G_REF2][i] == index[G_REF3][i] ||
 210                  index[G_REF3][i] == index[G_REF1][i] )
 211               {
 212                 fprintf(stderr,"Warning: reference triple (%d) will be skipped.\n",i);
 213                 fprintf(stderr,  "         index[1]: %d, index[2]: %d, index[3]: %d\n",
 214                         index[G_REF1][i],index[G_REF2][i],index[G_REF3][i]);
 215                 isize[NDX_REF1]--;
 216                 for (j=NDX_REF1; j<=NDX_REF3; j++)
 217                   srenew(index[j],isize[NDX_REF1]);
 218                 continue;
 219               }
 220             else /* everythings fine, copy that one */
 221               for (j=G_REF1; j<=G_REF3; j++)
 222                 index[j+4][i] = index[j][i];
 223         }
 224     }
 225   else if ( mode == 2 )
 226     {
 227       if ( isize[G_REF1] != isize[G_REF2] )
 228         gmx_fatal(FARGS,"For two particle SDF, reference groups 1 and 2"
 229                     "must have the same size.\n");
 230
 231
 232       for (i=0; i<isize[G_REF1]; i++)
 233         {
 234           /* check consistency for atoms 1 and 2 */
 235           for (j=G_REF1; j<=G_REF2; j++)
 236             ref_resind[j] = top.atoms.atom[index[j][i]].resind;
 237
 238
 239           if ( ref_resind[G_REF1] != ref_resind[G_REF2] ||
 240                index[G_REF1][i] == index[G_REF2][i] )
 241             {
 242               if ( ref_resind[G_REF1] != ref_resind[G_REF2] )
 243                 {
 244                   fprintf(stderr,"\nWarning: bond (%d) not from one molecule."
 245                           "Will not be used for SDF.\n",i);
 246                   fprintf(stderr,  "         resnr[1]: %d, resnr[2]: %d\n",
 247                           ref_resind[G_REF1],ref_resind[G_REF2]);
 248                 }
 249               else
 250                 {
 251                   fprintf(stderr,"\nWarning: atom1 and atom2 are identical."
 252                           "Bond (%d) will not be used for SDF.\n",i);
 253                   fprintf(stderr,  "         index[1]: %d, index[2]: %d\n",
 254                           index[G_REF1][i],index[G_REF2][i]);
 255                 }
 256               for (j=NDX_REF1; j<=NDX_REF2; j++)
 257                 {
 258                   for (k=i; k<isize[G_REF1]-1; k++)
 259                     index[j][k]=index[j][k+1];
 260                   isize[j]--;
 261                   srenew(index[j],isize[j]);
 262                 }
 263             }
 264         }
 265     }
 266
 267
 268   /* Read Atoms for refmol group */
 269   if ( bRef )
 270     {
 271       snew(index[G_REFMOL],1);
 272
 273
 274       for (i=G_REF1; i<=G_REF3; i++)
 275         ref_resind[i] = top.atoms.atom[index[i][0]].resind;
 276
 277
 278       for (i=0; i<natoms; i++)
 279         {
 280           if (  ref_resind[G_REF1] == top.atoms.atom[i].resind ||
 281                 ref_resind[G_REF2] == top.atoms.atom[i].resind ||
 282                 ref_resind[G_REF3] == top.atoms.atom[i].resind )
 283             nrefmol++;
 284         }
 285       srenew(index[G_REFMOL],nrefmol);
 286       isize[G_REFMOL] = nrefmol;
 287       nrefmol = 0;
 288
 289
 290       for (i=0; i<natoms; i++)
 291         {
 292           if (  ref_resind[G_REF1] == top.atoms.atom[i].resind ||
 293                 ref_resind[G_REF2] == top.atoms.atom[i].resind ||
 294                 ref_resind[G_REF3] == top.atoms.atom[i].resind )
 295             {
 296               index[G_REFMOL][nrefmol] = i;
 297               nrefmol++;
 298             }
 299         }
 300     }
 301
 302
 303   /* initialize some stuff */
 304   boxmin = min( norm(box[XX]), min( norm(box[YY]), norm(box[ZZ]) ) );
 305   hbox   = boxmin / 2.0;
 306
 307
 308   for (i=0; i<DIM; i++)
 309     {
 310       cutoff[i] = cutoff[i] / 2;
 311       nbin[i]   = (int)(2 * cutoff[i] / binwidth) + 1;
 312       invbinw = 1.0 / binwidth;
 313       tri_upper[i] = triangle[i] + dtri[i];
 314       tri_upper[i] = sqr(tri_upper[i]);
 315       tri_lower[i] = triangle[i] - dtri[i];
 316       tri_lower[i] = sqr(tri_lower[i]);
 317     }
 318
 319
 320   /* Allocate the array's for sdf */
 321   snew(count,nbin[0]+1);
 322   for(i=0; i<nbin[0]+1; i++)
 323     {
 324       snew(count[i],nbin[1]+1);
 325       for (j=0; j<nbin[1]+1; j++)
 326         snew(count[i][j],nbin[2]+1);
 327     }
 328
 329
 330   /* Allocate space for the coordinates */
 331   snew(x_i1,isize[G_SDF]);
 332   snew(x_refmol,isize[G_REFMOL]);
 333   for (i=0; i<isize[G_REFMOL]; i++)
 334     for (j=XX; j<=ZZ; j++)
 335       x_refmol[i][j] = 0;
 336
 337
 338   normfac = 0;
 339
 340
 341   do {
 342     /* Must init pbc every step because of pressure coupling */
 343     set_pbc(&pbc,ePBC,box);
 344     rm_pbc(&top.idef,ePBC,natoms,box,x,x);
 345
 346
 347     /* Dynamically build the ref triples */
 348     if ( mode == 2 )
 349       {
 350         isize[NDX_REF1]=0;
 351         for (j=NDX_REF1; j<=NDX_REF3; j++)
 352           srenew(index[j],isize[NDX_REF1]+1);
 353
 354
 355         /* consistancy of G_REF[1,2] has already been check */
 356         /* hence we can look for the third atom right away */
 357
 358
 359         for (i=0; i<isize[G_REF1]; i++)
 360           {
 361             for (j=0; j<isize[G_REF3]; j++)
 362               {
 363                 /* Avoid expensive stuff if possible */
 364                 if ( top.atoms.atom[index[G_REF1][i]].resind !=
 365                      top.atoms.atom[index[G_REF3][j]].resind &&
 366                      index[G_REF1][i] != index[G_REF3][j] &&
 367                      index[G_REF2][i] != index[G_REF3][j] )
 368                   {
 369                     pbc_dx(&pbc,x[index[G_REF1][i]],x[index[G_REF3][j]],dx);
 370                     delta = norm2(dx);
 371                     if ( delta < tri_upper[G_REF1] &&
 372                          delta > tri_lower[G_REF1] )
 373                       {
 374                         pbc_dx(&pbc,x[index[G_REF2][i]],x[index[G_REF3][j]],dx);
 375                         delta = norm2(dx);
 376                         if ( delta < tri_upper[G_REF2] &&
 377                              delta > tri_lower[G_REF2] )
 378                           {
 379                             /* found triple */
 380                             index[NDX_REF1][isize[NDX_REF1]]=index[G_REF1][i];
 381                             index[NDX_REF2][isize[NDX_REF1]]=index[G_REF2][i];
 382                             index[NDX_REF3][isize[NDX_REF1]]=index[G_REF3][j];
 383
 384
 385                             /* resize groups */
 386                             isize[NDX_REF1]++;
 387                             for (k=NDX_REF1; k<=NDX_REF3; k++)
 388                               srenew(index[k],isize[NDX_REF1]+1);
 389                           }
 390                       }
 391                   }
 392               }
 393           }
 394       }
 395     else if ( mode ==3 )
 396       {
 397         isize[NDX_REF1]=0;
 398         for (j=NDX_REF1; j<=NDX_REF3; j++)
 399           srenew(index[j],isize[NDX_REF1]+1);
 400
 401         /* consistancy will be checked while searching */
 402
 403
 404         for (i=0; i<isize[G_REF1]; i++)
 405           {
 406             for (j=0; j<isize[G_REF2]; j++)
 407               {
 408                 /* Avoid expensive stuff if possible */
 409                 if ( top.atoms.atom[index[G_REF1][i]].resind !=
 410                      top.atoms.atom[index[G_REF2][j]].resind &&
 411                      index[G_REF1][i] != index[G_REF2][j] )
 412                   {
 413                     pbc_dx(&pbc,x[index[G_REF1][i]],x[index[G_REF2][j]],dx);
 414                     delta = norm2(dx);
 415                     if ( delta < tri_upper[G_REF3] &&
 416                          delta > tri_lower[G_REF3] )
 417                       {
 418                         for (k=0; k<isize[G_REF3]; k++)
 419                           {
 420                             if ( top.atoms.atom[index[G_REF1][i]].resind !=
 421                                  top.atoms.atom[index[G_REF3][k]].resind &&
 422                                  top.atoms.atom[index[G_REF2][j]].resind !=
 423                                  top.atoms.atom[index[G_REF3][k]].resind &&
 424                                  index[G_REF1][i] != index[G_REF3][k] &&
 425                                  index[G_REF2][j] != index[G_REF3][k])
 426                               {
 427                                 pbc_dx(&pbc,x[index[G_REF1][i]],x[index[G_REF3][k]],dx);
 428                                 delta = norm2(dx);
 429                                 if ( delta < tri_upper[G_REF1] &&
 430                                      delta > tri_lower[G_REF1] )
 431                                   {
 432                                     pbc_dx(&pbc,x[index[G_REF2][j]],x[index[G_REF3][k]],dx);
 433                                     delta = norm2(dx);
 434                                     if ( delta < tri_upper[G_REF2] &&
 435                                          delta > tri_lower[G_REF2] )
 436                                       {
 437                                         /* found triple */
 438                                         index[NDX_REF1][isize[NDX_REF1]]=index[G_REF1][i];
 439                                         index[NDX_REF2][isize[NDX_REF1]]=index[G_REF2][j];
 440                                         index[NDX_REF3][isize[NDX_REF1]]=index[G_REF3][k];
 441
 442                                         /* resize groups */
 443                                         isize[NDX_REF1]++;
 444                                         for (l=NDX_REF1; l<=NDX_REF3; l++)
 445                                           srenew(index[l],isize[NDX_REF1]+1);
 446                                       }
 447                                   }
 448                               }
 449                           }
 450                       }
 451                   }
 452               }
 453           }
 454       }
 455
 456     for (i=0; i<isize[NDX_REF1]; i++)
 457       {
 458         /* setup the molecular coordinate system (i',j',k') */
 459         /* because the coodinate system of the box forms a unit matrix */
 460         /* (i',j',k') is identical with the rotation matrix */
 461         clear_mat(rot);
 462
 463
 464         /* k' = unitv(r(atom0) - r(atom1)) */
 465         pbc_dx(&pbc,x[index[NDX_REF1][i]],x[index[NDX_REF2][i]],k_mol);
 466         unitv(k_mol,rot[2]);
 467
 468         /* i' = unitv(k' x (r(atom2) - r(atom1))) */
 469         pbc_dx(&pbc,x[index[NDX_REF3][i]],x[index[NDX_REF2][i]],i1_mol);
 470         cprod(i1_mol,rot[2],i2_mol);
 471         unitv(i2_mol,rot[0]);
 472
 473         /* j' = k' x i' */
 474         cprod(rot[2],rot[0],rot[1]);
 475
 476
 477         /* set the point of reference */
 478         if ( mode == 2 )
 479           copy_rvec(x[index[NDX_REF3][i]],xi);
 480         else
 481           copy_rvec(x[index[NDX_REF1][i]],xi);
 482
 483
 484         /* make the reference */
 485         if ( bRef )
 486           {
 487             for (j=0; j<isize[G_REFMOL]; j++)
 488               {
 489                 pbc_dx(&pbc,xi,x[index[G_REFMOL][j]],dx);
 490                 mvmul(rot,dx,dx_mol);
 491                 rvec_inc(x_refmol[j],dx_mol);
 492                 for(k=XX; k<=ZZ; k++)
 493                    x_refmol[j][k] = x_refmol[j][k] / 2;
 494               }
 495           }
 496
 497
 498         /* Copy the indexed coordinates to a continuous array */
 499         for(j=0; j<isize[G_SDF]; j++)
 500           copy_rvec(x[index[G_SDF][j]],x_i1[j]);
 501
 502         /* count the SDF */
 503         for(j=0; j<isize[G_SDF]; j++)
 504           {
 505             /* Exclude peaks from the reference set */
 506             bInGroup=FALSE;
 507             for (k=NDX_REF1; k<=NDX_REF3; k++)
 508               if ( index[G_SDF][j] == index[k][i] )
 509                 bInGroup=TRUE;
 510
 511
 512             if ( !bInGroup )
 513               {
 514                 pbc_dx(&pbc,xi,x_i1[j],dx);
 515
 516                 /* transfer dx to the molecular coordinate system */
 517                 mvmul(rot,dx,dx_mol);
 518
 519
 520                 /* check cutoff's and count */
 521                 if ( dx_mol[XX] > -cutoff[XX] && dx_mol[XX] < cutoff[XX] )
 522                   if ( dx_mol[YY] > -cutoff[YY] && dx_mol[YY] < cutoff[YY] )
 523                     if ( dx_mol[ZZ] > -cutoff[ZZ] && dx_mol[ZZ] < cutoff[ZZ] )
 524                       {
 525                         X = (int)(floor(dx_mol[XX]*invbinw)) + (nbin[XX]-1)/2
 526 +1;
 527                         Y = (int)(floor(dx_mol[YY]*invbinw)) + (nbin[YY]-1)/2
 528 +1;
 529                         Z = (int)(floor(dx_mol[ZZ]*invbinw)) + (nbin[ZZ]-1)/2
 530 +1;
 531                         count[X][Y][Z]++;
 532                         normfac++;
 533                       }
 534               }
 535           }
 536       }
 537   } while (read_next_x(oenv,status,&t,natoms,x,box));
 538   fprintf(stderr,"\n");
 539
 540   close_trj(status);
 541
 542   sfree(x);
 543
 544
 545   /* write the reference strcture*/
 546   if ( bRef )
 547     {
 548       fp=ffopen(fnREF,"w");
 549       fprintf(fp,"%s\n",title);
 550       fprintf(fp,"  %d\n",isize[G_REFMOL]);
 551
 552
 553       for (i=0; i<isize[G_REFMOL]; i++)
 554         fprintf(fp,"%5d%5s%5s%5d%8.3f%8.3f%8.3f\n",
 555                 top.atoms.resinfo[top.atoms.atom[index[G_REFMOL][i]].resind].nr,
 556                 *(top.atoms.resinfo[top.atoms.atom[index[G_REFMOL][i]].resind].name),
 557                 *(top.atoms.atomname[index[G_REFMOL][i]]),i+1,
 558                 -1*x_refmol[i][XX],-1*x_refmol[i][YY],-1*x_refmol[i][ZZ]);
 559       /* Inserted -1* on the line above three times */
 560       fprintf(fp,"   10.00000   10.00000   10.00000\n");
 561       ffclose(fp);
 562       fprintf(stderr,"\nWrote reference structure. (%d Atoms)\n",isize[G_REFMOL]);
 563     }
 564
 565
 566   /* Calculate the mean probability density */
 567   fprintf(stderr,"\nNumber of configurations used for SDF: %d\n",(int)normfac);
 568
 569
 570   normfac = nbin[0]*nbin[1]*nbin[2] / normfac;
 571
 572
 573   fprintf(stderr,"\nMean probability density: %f\n",1/normfac);
 574
 575
 576   /* normalize the SDF and write output */
 577   /* see http://www.csc.fi/gopenmol/index.phtml for documentation */
 578   fp=ffopen(fnSDF,"wb");
 579
 580
 581   /* rank */
 582   i_write(fp,3);
 583
 584
 585   /* Type of surface */
 586   i_write(fp,42);
 587
 588
 589   /* Zdim, Ydim, Xdim */
 590   for (i=ZZ; i>=XX; i--)
 591     i_write(fp,nbin[i]);
 592
 593
 594   /* [Z,Y,X][min,max] (box corners in Angstroem)*/
 595   for (i=ZZ; i>=XX; i--)
 596     {
 597       f_write(fp,-cutoff[i]*10);
 598       f_write(fp,cutoff[i]*10);
 599     }
 600
 601
 602 /* Original Code
 603   for (i=1; i<nbin[2]+1; i++)
 604     for (j=1; j<nbin[1]+1; j++)
 605       for (k=1; k<nbin[0]+1; k++)
 606         {
 607           sdf = normfac * count[k][j][i];
 608           if ( sdf < min_sdf ) min_sdf = sdf;
 609           if ( sdf > max_sdf ) max_sdf = sdf;
 610           f_write(fp,sdf);
 611         }*/
 612 /* Changed Code to Mirror SDF to correct coordinates */
 613   for (i=nbin[2]; i>0; i--)
 614     for (j=nbin[1]; j>0; j--)
 615       for (k=nbin[0]; k>0; k--)
 616         {
 617           sdf = normfac * count[k][j][i];
 618           if ( sdf < min_sdf ) min_sdf = sdf;
 619           if ( sdf > max_sdf ) max_sdf = sdf;
 620           f_write(fp,sdf);
 621         }
 622
 623   fprintf(stderr,"\nMin: %f Max: %f\n",min_sdf,max_sdf);
 624
 625
 626   ffclose(fp);
 627
 628
 629   /* Give back the mem */
 630   for(i=0; i<nbin[0]+1; i++)
 631     {
 632       for (j=0; j<nbin[1]+1; j++)
 633         {
 634           sfree(count[i][j]);
 635         }
 636       sfree(count[i]);
 637     }
 638   sfree(count);
 639 }
 640
 641 int gmx_sdf(int argc,char *argv[])
 642 {
 643   const char *desc[] = {
 644     "g_sdf calculates the spatial distribution function (SDF) of a set of atoms",
 645     "within a coordinate system defined by three atoms. There is single body, ",
 646     "two body and three body SDF implemented (select with option -mode). ",
 647     "In the single body case the local coordinate system is defined by using",
 648     "a triple of atoms from one single molecule, for the two and three body case",
 649     "the configurations are dynamically searched complexes of two or three",
 650     "molecules (or residues) meeting certain distance consitions (see below).[PAR]",
 651     "The program needs a trajectory, a GROMACS run input file and an index ",
 652     "file to work. ",
 653     "You have to setup 4 groups in the index file before using g_sdf: [PAR]",
 654     "The first three groups are used to define the SDF coordinate system.",
 655     "The program will dynamically generate the atom triples according to ",
 656     "the selected -mode: ",
 657     "In -mode 1 the triples will be just the 1st, 2nd, 3rd, ... atoms from ",
 658     "groups 1, 2 and 3. Hence the nth entries in groups 1, 2 and 3 must be from the",
 659     "same residue. In -mode 2 the triples will be 1st, 2nd, 3rd, ... atoms from",
 660     "groups 1 and 2 (with the nth entries in groups 1 and 2 having the same res-id).",
 661     "For each pair from groups 1 and 2  group 3 is searched for an atom meeting the",
 662     "distance conditions set with -triangle and -dtri relative to atoms 1 and 2. In",
 663     "-mode 3 for each atom in group 1 group 2 is searched for an atom meeting the",
 664     "distance condition and if a pair is found group 3 is searched for an atom",
 665     "meeting the further conditions. The triple will only be used if all three atoms",
 666     "have different res-id's.[PAR]",
 667     "The local coordinate system is always defined using the following scheme:",
 668     "Atom 1 will be used as the point of origin for the SDF. ",
 669     "Atom 1 and 2 will define the principle axis (Z) of the coordinate system.",
 670     "The other two axis will be defined inplane (Y) and normal (X) to the plane through",
 671     "Atoms 1, 2 and 3. ",
 672     "The fourth group",
 673     "contains the atoms for which the SDF will be evaluated.[PAR]",
 674     "For -mode 2 and 3 you have to define the distance conditions for the ",
 675     "2 resp. 3 molecule complexes to be searched for using -triangle and -dtri.[PAR]",
 676     "The SDF will be sampled in cartesian coordinates.",
 677     "Use '-grid x y z' to define the size of the SDF grid around the ",
 678     "reference molecule. ",
 679     "The Volume of the SDF grid will be V=x*y*z (nm^3). ",
 680     "Use -bin to set the binwidth for grid.[PAR]",
 681     "The output will be a binary 3D-grid file (gom_plt.dat) in the .plt format that can be be",
 682     "read directly by gOpenMol. ",
 683     "The option -r will generate a .gro file with the reference molecule(s) transferred to",
 684     "the SDF coordinate system. Load this file into gOpenMol and display the",
 685     "SDF as a contour plot (see http://www.csc.fi/gopenmol/index.phtml for ",
 686     "further documentation). [PAR]",
 687     "For further information about SDF's have a look at: A. Vishnyakov, JPC A, 105,",
 688     "2001, 1702 and the references cited within."
 689   };
 690   output_env_t oenv;
 691   static bool bRef=FALSE;
 692   static int mode=1;
 693   static rvec triangle={0.0,0.0,0.0};
 694   static rvec dtri={0.03,0.03,0.03};
 695   static rvec cutoff={1,1,1};
 696   static real binwidth=0.05;
 697   t_pargs pa[] = {
 698     { "-mode",     FALSE, etINT, {&mode},
 699       "SDF in [1,2,3] particle mode"},
 700     { "-triangle", FALSE, etRVEC, {&triangle},
 701       "r(1,3), r(2,3), r(1,2)"},
 702     { "-dtri",     FALSE, etRVEC, {&dtri},
 703       "dr(1,3), dr(2,3), dr(1,2)"},
 704     { "-bin",      FALSE, etREAL, {&binwidth},
 705       "Binwidth for the 3D-grid (nm)" },
 706     { "-grid",      FALSE, etRVEC, {&cutoff},
 707       "Size of the 3D-grid (nm,nm,nm)"}
 708   };
 709 #define NPA asize(pa)
 710   const char       *fnTPS,*fnNDX,*fnREF;
 711
 712   t_filenm   fnm[] = {
 713     { efTRX, "-f",  NULL,     ffREAD },
 714     { efNDX, NULL,  NULL,     ffREAD },
 715     { efTPS, NULL,  NULL,     ffOPTRD },
 716     { efDAT, "-o",  "gom_plt",     ffWRITE },
 717     { efSTO, "-r",  "refmol",  ffOPTWR },
 718   };
 719 #define NFILE asize(fnm)
 720
 721   CopyRight(stderr,argv[0]);
 722   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_BE_NICE,
 723                     NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL,&oenv);
 724
 725
 726   fnTPS = ftp2fn_null(efTPS,NFILE,fnm);
 727   fnNDX = ftp2fn_null(efNDX,NFILE,fnm);
 728   fnREF = opt2fn_null("-r",NFILE,fnm);
 729   bRef  = opt2bSet("-r",NFILE,fnm);
 730
 731
 732
 733   if (!fnNDX)
 734     gmx_fatal(FARGS,"No index file specified\n"
 735                 "             Nothing to do!");
 736
 737
 738   if (!fnTPS)
 739     gmx_fatal(FARGS,"No topology file specified\n"
 740                 "             Nothing to do!");
 741
 742
 743   if ( bRef && (fn2ftp(fnREF) != efGRO))
 744     {
 745       fprintf(stderr,"\nOnly GROMACS format is supported for reference structures.\n");
 746       fprintf(stderr,"Option -r will be ignored!\n");
 747       bRef = FALSE;
 748     }
 749
 750
 751   if ((mode < 1) || (mode > 3))
 752     gmx_fatal(FARGS,"Wrong -mode selection. Chose 1-, 2- oder 3-particle mode.\n");
 753
 754
 755   do_sdf(fnNDX,fnTPS,ftp2fn(efTRX,NFILE,fnm),opt2fn("-o",NFILE,fnm),
 756          fnREF,bRef,cutoff,binwidth,mode,triangle,dtri,oenv);
 757
 758
 759   thanx(stderr);
 760
 761   return 0;
 762 }