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