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