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 "typedefs.h"
  24 #include "macros.h"
  25 #include "gromacs/math/vec.h"
  26 #include "pbc.h"
  27 #include "rmpbc.h"
  28 #include "copyrite.h"
  29 #include "gromacs/utility/futil.h"
  30 #include "gromacs/commandline/pargs.h"
  31 #include "gromacs/fileio/tpxio.h"
  32 #include "index.h"
  33 #include "gromacs/utility/smalloc.h"
  34 #include "nrnb.h"
  35 #include "gstat.h"
  36 #include "gromacs/fileio/matio.h"
  37 #include "gromacs/utility/fatalerror.h"
  38
  39
  40 #define G_REF1      0
  41 #define G_REF2      1
  42 #define G_REF3      2
  43 #define G_SDF       3
  44 #define NDX_REF1    4
  45 #define NDX_REF2    5
  46 #define NDX_REF3    6
  47 #define G_REFMOL    7
  48
  49 static void i_write(FILE *output, int value)
  50 {
  51   if(fwrite(&value,sizeof(int),1,output) != 1)
  52   {
  53     gmx_fatal(FARGS,"Error writing to output file");
  54   }
  55 }
  56
  57
  58 static void f_write(FILE *output,float value)
  59 {
  60   if(fwrite(&value,sizeof(float),1,output) != 1)
  61   {
  62     gmx_fatal(FARGS,"Error writing to output file");
  63   }
  64 }
  65
  66
  67 static void do_sdf(const char *fnNDX,const char *fnTPS,const char *fnTRX,
  68                    const char *fnSDF, const char *fnREF, gmx_bool bRef,
  69                    rvec cutoff, real binwidth, int mode, rvec triangle,
  70                    rvec dtri, const output_env_t oenv)
  71 {
  72   FILE       *fp;
  73   t_trxstatus *status;
  74   int        ng,natoms,i,j,k,l,X,Y,Z,lc,dest;
  75   ivec       nbin;
  76   int        ***count;
  77   /* real       ***sdf; */
  78   real       sdf,min_sdf=1e10,max_sdf=0;
  79   char       **grpname;
  80   int        *isize;
  81   int        isize_cg=0;
  82   int        isize_ref=3;
  83   int        ref_resind[3]={0};
  84   int        nrefmol=0,refc=0;
  85   atom_id    **index;
  86   atom_id    *index_cg=NULL;
  87   atom_id    *index_ref=NULL;
  88   real       t,boxmin,hbox,normfac;
  89   real       invbinw;
  90   rvec       tri_upper,tri_lower;
  91   rvec       *x,xcog,dx,*x_i1,xi,*x_refmol;
  92   matrix     box;
  93   matrix     rot; /* rotation matrix := unit vectors for the molecule frame */
  94   rvec       k_mol,i1_mol,i2_mol,dx_mol;
  95   real       delta;
  96   atom_id    ix,jx;
  97   t_topology top;
  98   gmx_rmpbc_t  gpbc=NULL;
  99   int        ePBC=-1;
 100   t_pbc      pbc;
 101   gmx_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   gpbc = gmx_rmpbc_init(&top.idef,ePBC,natoms,box);
 341
 342   do {
 343     /* Must init pbc every step because of pressure coupling */
 344     set_pbc(&pbc,ePBC,box);
 345     gmx_rmpbc(gpbc,natoms,box,x);
 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   gmx_rmpbc_done(gpbc);
 541
 542
 543   close_trj(status);
 544
 545   sfree(x);
 546
 547
 548   /* write the reference strcture*/
 549   if ( bRef )
 550     {
 551       fp=gmx_ffopen(fnREF,"w");
 552       fprintf(fp,"%s\n",title);
 553       fprintf(fp,"  %d\n",isize[G_REFMOL]);
 554
 555
 556       for (i=0; i<isize[G_REFMOL]; i++)
 557         fprintf(fp,"%5d%5s%5s%5d%8.3f%8.3f%8.3f\n",
 558                 top.atoms.resinfo[top.atoms.atom[index[G_REFMOL][i]].resind].nr,
 559                 *(top.atoms.resinfo[top.atoms.atom[index[G_REFMOL][i]].resind].name),
 560                 *(top.atoms.atomname[index[G_REFMOL][i]]),i+1,
 561                 -1*x_refmol[i][XX],-1*x_refmol[i][YY],-1*x_refmol[i][ZZ]);
 562       /* Inserted -1* on the line above three times */
 563       fprintf(fp,"   10.00000   10.00000   10.00000\n");
 564       gmx_ffclose(fp);
 565       fprintf(stderr,"\nWrote reference structure. (%d Atoms)\n",isize[G_REFMOL]);
 566     }
 567
 568
 569   /* Calculate the mean probability density */
 570   fprintf(stderr,"\nNumber of configurations used for SDF: %d\n",(int)normfac);
 571
 572
 573   normfac = nbin[0]*nbin[1]*nbin[2] / normfac;
 574
 575
 576   fprintf(stderr,"\nMean probability density: %f\n",1/normfac);
 577
 578
 579   /* normalize the SDF and write output */
 580   /* see http://www.csc.fi/gopenmol/index.phtml for documentation */
 581   fp=gmx_ffopen(fnSDF,"wb");
 582
 583
 584   /* rank */
 585   i_write(fp,3);
 586
 587
 588   /* Type of surface */
 589   i_write(fp,42);
 590
 591
 592   /* Zdim, Ydim, Xdim */
 593   for (i=ZZ; i>=XX; i--)
 594     i_write(fp,nbin[i]);
 595
 596
 597   /* [Z,Y,X][min,max] (box corners in Angstroem)*/
 598   for (i=ZZ; i>=XX; i--)
 599     {
 600       f_write(fp,-cutoff[i]*10);
 601       f_write(fp,cutoff[i]*10);
 602     }
 603
 604
 605 /* Original Code
 606   for (i=1; i<nbin[2]+1; i++)
 607     for (j=1; j<nbin[1]+1; j++)
 608       for (k=1; k<nbin[0]+1; k++)
 609         {
 610           sdf = normfac * count[k][j][i];
 611           if ( sdf < min_sdf ) min_sdf = sdf;
 612           if ( sdf > max_sdf ) max_sdf = sdf;
 613           f_write(fp,sdf);
 614         }*/
 615 /* Changed Code to Mirror SDF to correct coordinates */
 616   for (i=nbin[2]; i>0; i--)
 617     for (j=nbin[1]; j>0; j--)
 618       for (k=nbin[0]; k>0; k--)
 619         {
 620           sdf = normfac * count[k][j][i];
 621           if ( sdf < min_sdf ) min_sdf = sdf;
 622           if ( sdf > max_sdf ) max_sdf = sdf;
 623           f_write(fp,sdf);
 624         }
 625
 626   fprintf(stderr,"\nMin: %f Max: %f\n",min_sdf,max_sdf);
 627
 628
 629   gmx_ffclose(fp);
 630
 631
 632   /* Give back the mem */
 633   for(i=0; i<nbin[0]+1; i++)
 634     {
 635       for (j=0; j<nbin[1]+1; j++)
 636         {
 637           sfree(count[i][j]);
 638         }
 639       sfree(count[i]);
 640     }
 641   sfree(count);
 642 }
 643
 644 int gmx_sdf(int argc,char *argv[])
 645 {
 646   const char *desc[] = {
 647     "[TT]g_sdf[tt] calculates the spatial distribution function (SDF) of a set of atoms",
 648     "within a coordinate system defined by three atoms. There is single body, ",
 649     "two body and three body SDF implemented (select with option [TT]-mode[tt]). ",
 650     "In the single body case the local coordinate system is defined by using",
 651     "a triple of atoms from one single molecule, for the two and three body case",
 652     "the configurations are dynamically searched complexes of two or three",
 653     "molecules (or residues) meeting certain distance consitions (see below).[PAR]",
 654     "The program needs a trajectory, a GROMACS run input file and an index ",
 655     "file to work. ",
 656     "You have to setup 4 groups in the index file before using g_sdf: [PAR]",
 657     "The first three groups are used to define the SDF coordinate system.",
 658     "The program will dynamically generate the atom triples according to ",
 659     "the selected [TT]-mode[tt]: ",
 660     "In [TT]-mode[tt] 1 the triples will be just the 1st, 2nd, 3rd, ... atoms from ",
 661     "groups 1, 2 and 3. Hence the nth entries in groups 1, 2 and 3 must be from the",
 662     "same residue. In [TT]-mode[tt] 2 the triples will be 1st, 2nd, 3rd, ... atoms from",
 663     "groups 1 and 2 (with the nth entries in groups 1 and 2 having the same res-id).",
 664     "For each pair from groups 1 and 2  group 3 is searched for an atom meeting the",
 665     "distance conditions set with [TT]-triangle[tt] and [TT]-dtri[tt] relative to atoms 1 and 2. In",
 666     "[TT]-mode[tt] 3 for each atom in group 1 group 2 is searched for an atom meeting the",
 667     "distance condition and if a pair is found group 3 is searched for an atom",
 668     "meeting the further conditions. The triple will only be used if all three atoms",
 669     "have different res-id's.[PAR]",
 670     "The local coordinate system is always defined using the following scheme:",
 671     "Atom 1 will be used as the point of origin for the SDF. ",
 672     "Atom 1 and 2 will define the principle axis (Z) of the coordinate system.",
 673     "The other two axis will be defined inplane (Y) and normal (X) to the plane through",
 674     "Atoms 1, 2 and 3. ",
 675     "The fourth group",
 676     "contains the atoms for which the SDF will be evaluated.[PAR]",
 677     "For [TT]-mode[tt] 2 and 3 you have to define the distance conditions for the ",
 678     "2 resp. 3 molecule complexes to be searched for using [TT]-triangle[tt] and [TT]-dtri[tt].[PAR]",
 679     "The SDF will be sampled in cartesian coordinates.",
 680     "Use [TT]-grid x y z[tt] to define the size of the SDF grid around the ",
 681     "reference molecule. ",
 682     "The Volume of the SDF grid will be V=x*y*z (nm^3). ",
 683     "Use [TT]-bin[tt] to set the binwidth for grid.[PAR]",
 684     "The output will be a binary 3D-grid file ([TT]gom_plt.dat[tt]) in the [TT].plt[tt] format that can be be",
 685     "read directly by gOpenMol. ",
 686     "The option [TT]-r[tt] will generate a [TT].gro[tt] file with the reference molecule(s) transferred to",
 687     "the SDF coordinate system. Load this file into gOpenMol and display the",
 688     "SDF as a contour plot (see http://www.csc.fi/gopenmol/index.phtml for ",
 689     "further documentation). [PAR]",
 690     "For further information about SDF's have a look at: A. Vishnyakov, JPC A, 105,",
 691     "2001, 1702 and the references cited within."
 692   };
 693   output_env_t oenv;
 694   static gmx_bool bRef=FALSE;
 695   static int mode=1;
 696   static rvec triangle={0.0,0.0,0.0};
 697   static rvec dtri={0.03,0.03,0.03};
 698   static rvec cutoff={1,1,1};
 699   static real binwidth=0.05;
 700   t_pargs pa[] = {
 701     { "-mode",     FALSE, etINT, {&mode},
 702       "SDF in [1,2,3] particle mode"},
 703     { "-triangle", FALSE, etRVEC, {&triangle},
 704       "r(1,3), r(2,3), r(1,2)"},
 705     { "-dtri",     FALSE, etRVEC, {&dtri},
 706       "dr(1,3), dr(2,3), dr(1,2)"},
 707     { "-bin",      FALSE, etREAL, {&binwidth},
 708       "Binwidth for the 3D-grid (nm)" },
 709     { "-grid",      FALSE, etRVEC, {&cutoff},
 710       "Size of the 3D-grid (nm,nm,nm)"}
 711   };
 712 #define NPA asize(pa)
 713   const char       *fnTPS,*fnNDX,*fnREF;
 714
 715   t_filenm   fnm[] = {
 716     { efTRX, "-f",  NULL,     ffREAD },
 717     { efNDX, NULL,  NULL,     ffREAD },
 718     { efTPS, NULL,  NULL,     ffOPTRD },
 719     { efDAT, "-o",  "gom_plt",     ffWRITE },
 720     { efSTO, "-r",  "refmol",  ffOPTWR },
 721   };
 722 #define NFILE asize(fnm)
 723
 724   CopyRight(stderr,argv[0]);
 725   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_BE_NICE,
 726                     NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL,&oenv);
 727
 728
 729   fnTPS = ftp2fn_null(efTPS,NFILE,fnm);
 730   fnNDX = ftp2fn_null(efNDX,NFILE,fnm);
 731   fnREF = opt2fn_null("-r",NFILE,fnm);
 732   bRef  = opt2bSet("-r",NFILE,fnm);
 733
 734
 735
 736   if (!fnNDX)
 737     gmx_fatal(FARGS,"No index file specified\n"
 738                 "             Nothing to do!");
 739
 740
 741   if (!fnTPS)
 742     gmx_fatal(FARGS,"No topology file specified\n"
 743                 "             Nothing to do!");
 744
 745
 746   if ( bRef && (fn2ftp(fnREF) != efGRO))
 747     {
 748       fprintf(stderr,"\nOnly GROMACS format is supported for reference structures.\n");
 749       fprintf(stderr,"Option -r will be ignored!\n");
 750       bRef = FALSE;
 751     }
 752
 753
 754   if ((mode < 1) || (mode > 3))
 755     gmx_fatal(FARGS,"Wrong -mode selection. Chose 1-, 2- oder 3-particle mode.\n");
 756
 757
 758   do_sdf(fnNDX,fnTPS,ftp2fn(efTRX,NFILE,fnm),opt2fn("-o",NFILE,fnm),
 759          fnREF,bRef,cutoff,binwidth,mode,triangle,dtri,oenv);
 760
 761
 762   gmx_thanx(stderr);
 763
 764   return 0;
 765 }
 766
 767 int
 768 main(int argc, char *argv[])
 769 {
 770   gmx_sdf(argc,argv);
 771   return 0;
 772 }