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