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