src/tools/gmx_dyecoupl.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  10  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  11  * Copyright (c) 2001-2009, The GROMACS development team,
  12  * check out http://www.gromacs.org for more information.
  13
  14  * This program is free software; you can redistribute it and/or
  15  * modify it under the terms of the GNU General Public License
  16  * as published by the Free Software Foundation; either version 2
  17  * of the License, or (at your option) any later version.
  18  *
  19  * If you want to redistribute modifications, please consider that
  20  * scientific software is very special. Version control is crucial -
  21  * bugs must be traceable. We will be happy to consider code for
  22  * inclusion in the official distribution, but derived work must not
  23  * be called official GROMACS. Details are found in the README & COPYING
  24  * files - if they are missing, get the official version at www.gromacs.org.
  25  *
  26  * To help us fund GROMACS development, we humbly ask that you cite
  27  * the papers on the package - you can find them in the top README file.
  28  *
  29  * For more info, check our website at http://www.gromacs.org
  30  */
  31 #include <copyrite.h>
  32 #include <filenm.h>
  33 #include <macros.h>
  34 #include <pbc.h>
  35 #include <smalloc.h>
  36 #include <statutil.h>
  37 #include <vec.h>
  38 #include <xvgr.h>
  39
  40
  41 int gmx_dyecoupl(int argc, char *argv[])
  42 {
  43     const char *desc[] =
  44     {
  45             "This tool extracts dye dynamics from trajectory files.",
  46             "Currently, R and kappa^2 between dyes is extracted for (F)RET",
  47             "simulations with assumed dipolar coupling as in the Foerster equation.",
  48             "It further allows the calculation of R(t) and kappa^2(t), R and",
  49             "kappa^2 histograms and averages, as well as the instantaneous FRET",
  50             "efficiency E(t) for a specified Foerster radius R_0 (switch [TT]-R0[tt]).",
  51             "The input dyes have to be whole (see res and mol pbc options",
  52             "in [TT]trjconv[tt]).",
  53             "The dye transition dipole moment has to be defined by at least",
  54             "a single atom pair, however multiple atom pairs can be provided ",
  55             "in the index file. The distance R is calculated on the basis of",
  56             "the COMs of the given atom pairs.",
  57             "The [TT]-pbcdist[tt] option calculates distances to the nearest periodic",
  58             "image instead to the distance in the box. This works however only,"
  59             "for periodic boundaries in all 3 dimensions.",
  60             "The [TT]-norm[tt] option (area-) normalizes the histograms."
  61     };
  62
  63         static gmx_bool bPBCdist = FALSE, bNormHist = FALSE;
  64     int histbins = 50;
  65     output_env_t oenv;
  66     real R0=-1;
  67
  68     t_pargs pa[] =
  69     {
  70             { "-pbcdist", FALSE, etBOOL, { &bPBCdist },"Distance R based on PBC" },
  71             { "-norm", FALSE, etBOOL, { &bNormHist },"Normalize histograms" },
  72             { "-bins", FALSE, etINT, {&histbins},"# of histogram bins" },
  73             { "-R0", FALSE, etREAL, {&R0},"Foerster radius including kappa^2=2/3 in nm" }
  74     };
  75 #define NPA asize(pa)
  76
  77     t_filenm fnm[] =
  78     {
  79             { efTRX, "-f", NULL, ffREAD },
  80             { efNDX, NULL, NULL, ffREAD },
  81             { efXVG, "-ot", "rkappa",ffOPTWR },
  82             { efXVG, "-oe", "insteff",ffOPTWR },
  83             { efDAT, "-o", "rkappa",ffOPTWR },
  84             { efXVG, "-rhist","rhist", ffOPTWR },
  85             { efXVG, "-khist", "khist", ffOPTWR }
  86     };
  87 #define NFILE asize(fnm)
  88
  89
  90     const char *in_trajfile, *in_ndxfile, *out_xvgrkfile = NULL, *out_xvginstefffile = NULL, *out_xvgrhistfile = NULL, *out_xvgkhistfile = NULL,*out_datfile=NULL;
  91     gmx_bool bHaveFirstFrame, bHaveNextFrame, indexOK = TRUE;
  92     int ndon, nacc;
  93     atom_id *donindex, *accindex;
  94     char *grpnm;
  95     t_atoms *atoms = NULL;
  96     t_trxstatus *status;
  97     t_trxframe fr;
  98
  99     int flags;
 100     int allocblock = 1000;
 101     real histexpand = 1e-6;
 102     rvec donvec, accvec, donpos, accpos, dist, distnorm;
 103     int natoms;
 104
 105     /*we rely on PBC autodetection (...currently)*/
 106     int ePBC = -1;
 107
 108     real *rvalues=NULL, *kappa2values=NULL, *rhist=NULL, *khist=NULL;
 109     t_pbc *pbc=NULL;
 110     int i, bin;
 111     FILE *rkfp = NULL, *rhfp = NULL, *khfp = NULL,*datfp=NULL,*iefp=NULL;
 112     gmx_bool bRKout, bRhistout, bKhistout,bDatout,bInstEffout,grident;
 113
 114     const char *rkleg[2] = { "R", "\\f{Symbol}k\\f{}\\S2\\N" };
 115     const char *rhleg[1] = { "p(R)" };
 116     const char *khleg[1] = { "p(\\f{Symbol}k\\f{}\\S2\\N)" };
 117     const char *ieleg[1] = { "E\\sRET\\N(t)" };
 118
 119     real R, kappa2, insteff, Rs = 0., kappa2s = 0., insteffs=0., rmax, rmin, kmin = 0., kmax = 4.,
 120             rrange, krange, rincr, kincr,Rfrac;
 121     int rkcount = 0, rblocksallocated = 0, kblocksallocated = 0;
 122
 123     CopyRight(stderr, argv[0]);
 124     parse_common_args(&argc,argv,PCA_CAN_BEGIN | PCA_CAN_END | PCA_CAN_VIEW | PCA_TIME_UNIT | PCA_BE_NICE, NFILE,fnm,NPA,pa,asize(desc),desc, 0,NULL,&oenv);
 125
 126
 127     /* Check command line options for filenames and set bool flags when switch used*/
 128     in_trajfile = opt2fn("-f", NFILE, fnm);
 129     in_ndxfile = opt2fn("-n", NFILE, fnm);
 130     out_xvgrkfile = opt2fn("-ot", NFILE, fnm);
 131     out_xvgrhistfile = opt2fn("-rhist", NFILE, fnm);
 132     out_xvgkhistfile = opt2fn("-khist", NFILE, fnm);
 133     out_xvginstefffile = opt2fn("-oe", NFILE, fnm);
 134     out_datfile = opt2fn("-o",NFILE,fnm);
 135
 136     bRKout = opt2bSet("-ot", NFILE, fnm);
 137     bRhistout = opt2bSet("-rhist", NFILE, fnm);
 138     bKhistout = opt2bSet("-khist", NFILE, fnm);
 139     bDatout = opt2bSet("-o", NFILE, fnm);
 140     bInstEffout = opt2bSet("-oe", NFILE, fnm);
 141
 142
 143     /* PBC warning. */
 144     if (bPBCdist)
 145     {
 146         printf("Calculating distances to periodic image.\n");
 147         printf("Be careful! This produces only valid results for PBC in all three dimensions\n");
 148     }
 149
 150
 151     if (bInstEffout && R0<=0.)
 152     {
 153         gmx_fatal(FARGS,"You have to specify R0 and R0 has to be larger than 0 nm.\n\n");
 154     }
 155
 156     printf("Select group with donor atom pairs defining the transition moment\n");
 157     get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &ndon, &donindex,&grpnm);
 158
 159     printf("Select group with acceptor atom pairs defining the transition moment\n");
 160     get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &nacc, &accindex,&grpnm);
 161
 162     /*check if groups are identical*/
 163     grident=TRUE;
 164
 165     if (ndon == nacc)
 166     {
 167         for (i=0;i<nacc;i++)
 168         {
 169             if(accindex[i] != donindex[i])
 170             {
 171                 grident=FALSE;
 172                 break;
 173             }
 174         }
 175     }
 176
 177     if (grident)
 178     {
 179         gmx_fatal(FARGS,"Donor and acceptor group are identical. This makes no sense.");
 180     }
 181
 182     printf("Reading first frame\n");
 183     /* open trx file for reading */
 184     flags=0;
 185     flags = flags | TRX_READ_X;
 186     bHaveFirstFrame = read_first_frame(oenv, &status, in_trajfile, &fr, flags);
 187
 188     if (bHaveFirstFrame)
 189     {
 190         printf("First frame is OK\n");
 191         natoms = fr.natoms;
 192         if ((ndon % 2 != 0) || (nacc % 2 != 0))
 193         {
 194             indexOK = FALSE;
 195         }
 196         else
 197         {
 198             for (i = 0; i < ndon;i++)
 199             {
 200                 if (donindex[i] >= natoms)
 201                     indexOK = FALSE;
 202             }
 203             for (i = 0; i < nacc;i++)
 204             {
 205                 if (accindex[i] >= natoms)
 206                     indexOK = FALSE;
 207             }
 208         }
 209
 210         if (indexOK)
 211         {
 212
 213             if (bDatout)
 214             {
 215                 datfp = fopen(out_datfile,"w");
 216             }
 217
 218             if (bRKout)
 219             {
 220                 rkfp = xvgropen(out_xvgrkfile,
 221                         "Distance and \\f{Symbol}k\\f{}\\S2\\N trajectory",
 222                         "Time (ps)", "Distance (nm) / \\f{Symbol}k\\f{}\\S2\\N",
 223                         oenv);
 224                 xvgr_legend(rkfp, 2, rkleg, oenv);
 225             }
 226
 227             if (bInstEffout)
 228             {
 229                 iefp = xvgropen(out_xvginstefffile,
 230                         "Instantaneous RET Efficiency",
 231                         "Time (ps)", "RET Efficiency",
 232                         oenv);
 233                 xvgr_legend(iefp, 1, ieleg, oenv);
 234             }
 235
 236
 237             if (bRhistout)
 238             {
 239                 snew(rvalues, allocblock);
 240                 rblocksallocated += 1;
 241                 snew(rhist, histbins);
 242             }
 243
 244             if (bKhistout)
 245             {
 246                 snew(kappa2values, allocblock);
 247                 kblocksallocated += 1;
 248                 snew(khist, histbins);
 249             }
 250
 251             do
 252             {
 253                 clear_rvec(donvec);
 254                 clear_rvec(accvec);
 255                 clear_rvec(donpos);
 256                 clear_rvec(accpos);
 257                 for (i = 0; i < ndon / 2; i++)
 258                 {
 259                     rvec_sub(donvec, fr.x[donindex[2 * i]], donvec);
 260                     rvec_add(donvec, fr.x[donindex[2 * i + 1]], donvec);
 261                     rvec_add(donpos, fr.x[donindex[2 * i]], donpos);
 262                     rvec_add(donpos, fr.x[donindex[2 * i + 1]], donpos);
 263                 }
 264
 265                 for (i = 0; i < nacc / 2; i++)
 266                 {
 267                     rvec_sub(accvec, fr.x[accindex[2 * i]], accvec);
 268                     rvec_add(accvec, fr.x[accindex[2 * i + 1]], accvec);
 269                     rvec_add(accpos, fr.x[accindex[2 * i]], accpos);
 270                     rvec_add(accpos, fr.x[accindex[2 * i + 1]], accpos);
 271                 }
 272
 273                 unitv(donvec, donvec);
 274                 unitv(accvec, accvec);
 275
 276                 svmul((real) 1. / ndon, donpos, donpos);
 277                 svmul((real) 1. / nacc, accpos, accpos);
 278
 279                 if (bPBCdist)
 280                 {
 281                     set_pbc(pbc, ePBC, fr.box);
 282                     pbc_dx(pbc, donpos, accpos, dist);
 283                 }
 284                 else
 285                 {
 286                     rvec_sub(donpos, accpos, dist);
 287                 }
 288
 289                 unitv(dist, distnorm);
 290                 R = norm(dist);
 291                 kappa2 = iprod(donvec, accvec)- 3.* (iprod(donvec, distnorm) * iprod(distnorm, accvec));
 292                 kappa2 *= kappa2;
 293                 if (R0>0)
 294                 {
 295                     Rfrac=R/R0;
 296                     insteff=1/(1+(Rfrac*Rfrac*Rfrac*Rfrac*Rfrac*Rfrac)*2/3/kappa2);
 297                     insteffs+=insteff;
 298
 299                     if (bInstEffout)
 300                     {
 301                         fprintf(iefp, "%12.7f %12.7f\n", fr.time, insteff);
 302                     }
 303                 }
 304
 305
 306                 Rs += R;
 307                 kappa2s += kappa2;
 308                 rkcount++;
 309
 310                 if (bRKout)
 311                     fprintf(rkfp, "%12.7f %12.7f %12.7f\n", fr.time, R, kappa2);
 312
 313                 if (bDatout)
 314                     fprintf(datfp, "%12.7f %12.7f %12.7f\n", fr.time, R, kappa2);
 315
 316                 if (bRhistout)
 317                 {
 318                     rvalues[rkcount-1] = R;
 319                     if (rkcount % allocblock == 0)
 320                     {
 321                         srenew(rvalues, allocblock*(rblocksallocated+1));
 322                         rblocksallocated += 1;
 323                     }
 324                 }
 325
 326                 if (bKhistout)
 327                 {
 328                     kappa2values[rkcount-1] = kappa2;
 329                     if (rkcount % allocblock == 0)
 330                     {
 331                         srenew(kappa2values, allocblock*(kblocksallocated+1));
 332                         kblocksallocated += 1;
 333                     }
 334                 }
 335
 336                 bHaveNextFrame = read_next_frame(oenv, status, &fr);
 337             } while (bHaveNextFrame);
 338
 339             if (bRKout)
 340                 ffclose(rkfp);
 341
 342             if (bDatout)
 343                 ffclose(datfp);
 344
 345             if (bInstEffout)
 346                 ffclose(iefp);
 347
 348
 349             if (bRhistout)
 350             {
 351                 printf("Writing R-Histogram\n");
 352                 rmin = rvalues[0];
 353                 rmax = rvalues[0];
 354                 for (i = 1; i < rkcount; i++)
 355                 {
 356                     if (rvalues[i] < rmin)
 357                         rmin = rvalues[i];
 358                     else if (rvalues[i] > rmax)
 359                         rmax = rvalues[i];
 360                 }
 361                 rmin -= histexpand;
 362                 rmax += histexpand;
 363
 364                 rrange = rmax - rmin;
 365                 rincr = rrange / histbins;
 366
 367                 for (i = 1; i < rkcount; i++)
 368                 {
 369                     bin = (int) ((rvalues[i] - rmin) / rincr);
 370                     rhist[bin] += 1;
 371                 }
 372                 if (bNormHist)
 373                 {
 374                     for (i = 0; i < histbins; i++)
 375                         rhist[i] /= rkcount * rrange/histbins;
 376                     rhfp = xvgropen(out_xvgrhistfile, "Distance Distribution",
 377                             "R (nm)", "Normalized Probability", oenv);
 378                 } else
 379                 {
 380                     rhfp = xvgropen(out_xvgrhistfile, "Distance Distribution",
 381                             "R (nm)", "Probability", oenv);
 382                 }
 383                 xvgr_legend(rhfp, 1, rhleg, oenv);
 384                 for (i = 0; i < histbins; i++)
 385                 {
 386                     fprintf(rhfp, "%12.7f %12.7f\n", (i + 0.5) * rincr + rmin,
 387                             rhist[i]);
 388                 }
 389                 ffclose(rhfp);
 390             }
 391
 392             if (bKhistout)
 393             {
 394                 printf("Writing kappa^2-Histogram\n");
 395                 krange = kmax - kmin;
 396                 kincr = krange / histbins;
 397
 398                 for (i = 1; i < rkcount; i++)
 399                 {
 400                     bin = (int) ((kappa2values[i] - kmin) / kincr);
 401                     khist[bin] += 1;
 402                 }
 403                 if (bNormHist)
 404                 {
 405                     for (i = 0; i < histbins; i++)
 406                         khist[i] /= rkcount * krange/histbins;
 407                     khfp = xvgropen(out_xvgkhistfile,
 408                             "\\f{Symbol}k\\f{}\\S2\\N Distribution",
 409                             "\\f{Symbol}k\\f{}\\S2\\N",
 410                             "Normalized Probability", oenv);
 411                 } else
 412                 {
 413                     khfp = xvgropen(out_xvgkhistfile,
 414                             "\\f{Symbol}k\\f{}\\S2\\N Distribution",
 415                             "\\f{Symbol}k\\f{}\\S2\\N", "Probability", oenv);
 416                 }
 417                 xvgr_legend(khfp, 1, khleg, oenv);
 418                 for (i = 0; i < histbins; i++)
 419                 {
 420                     fprintf(khfp, "%12.7f %12.7f\n", (i + 0.5) * kincr + kmin,
 421                             khist[i]);
 422                 }
 423                 ffclose(khfp);
 424             }
 425
 426             printf("\nAverages:\n");
 427             printf("R_avg   = %8.4f nm\nKappa^2 = %8.4f\n", Rs / rkcount,
 428                     kappa2s / rkcount);
 429             if (R0>0)
 430             {
 431                 printf("E_RETavg   = %8.4f\n", insteffs / rkcount);
 432             }
 433             please_cite(stdout,"Hoefling2011");
 434         }
 435         else
 436         {
 437             gmx_fatal(FARGS,"Index file invalid, check your index file for correct pairs.\n");
 438         }
 439     }
 440     else
 441     {
 442         gmx_fatal(FARGS,"Could not read first frame of the trajectory.\n");
 443     }
 444
 445     thanx(stderr);
 446     return 0;
 447 }
 448