src/gromacs/gmxana/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     if (!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))
 124     {
 125         return 0;
 126     }
 127
 128
 129     /* Check command line options for filenames and set bool flags when switch used*/
 130     in_trajfile        = opt2fn("-f", NFILE, fnm);
 131     in_ndxfile         = opt2fn("-n", NFILE, fnm);
 132     out_xvgrkfile      = opt2fn("-ot", NFILE, fnm);
 133     out_xvgrhistfile   = opt2fn("-rhist", NFILE, fnm);
 134     out_xvgkhistfile   = opt2fn("-khist", NFILE, fnm);
 135     out_xvginstefffile = opt2fn("-oe", NFILE, fnm);
 136     out_datfile        = opt2fn("-o", NFILE, fnm);
 137
 138     bRKout      = opt2bSet("-ot", NFILE, fnm);
 139     bRhistout   = opt2bSet("-rhist", NFILE, fnm);
 140     bKhistout   = opt2bSet("-khist", NFILE, fnm);
 141     bDatout     = opt2bSet("-o", NFILE, fnm);
 142     bInstEffout = opt2bSet("-oe", NFILE, fnm);
 143
 144
 145     /* PBC warning. */
 146     if (bPBCdist)
 147     {
 148         printf("Calculating distances to periodic image.\n");
 149         printf("Be careful! This produces only valid results for PBC in all three dimensions\n");
 150     }
 151
 152
 153     if (bInstEffout && R0 <= 0.)
 154     {
 155         gmx_fatal(FARGS, "You have to specify R0 and R0 has to be larger than 0 nm.\n\n");
 156     }
 157
 158     printf("Select group with donor atom pairs defining the transition moment\n");
 159     get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &ndon, &donindex, &grpnm);
 160
 161     printf("Select group with acceptor atom pairs defining the transition moment\n");
 162     get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &nacc, &accindex, &grpnm);
 163
 164     /*check if groups are identical*/
 165     grident = TRUE;
 166
 167     if (ndon == nacc)
 168     {
 169         for (i = 0; i < nacc; i++)
 170         {
 171             if (accindex[i] != donindex[i])
 172             {
 173                 grident = FALSE;
 174                 break;
 175             }
 176         }
 177     }
 178
 179     if (grident)
 180     {
 181         gmx_fatal(FARGS, "Donor and acceptor group are identical. This makes no sense.");
 182     }
 183
 184     printf("Reading first frame\n");
 185     /* open trx file for reading */
 186     flags           = 0;
 187     flags           = flags | TRX_READ_X;
 188     bHaveFirstFrame = read_first_frame(oenv, &status, in_trajfile, &fr, flags);
 189
 190     if (bHaveFirstFrame)
 191     {
 192         printf("First frame is OK\n");
 193         natoms = fr.natoms;
 194         if ((ndon % 2 != 0) || (nacc % 2 != 0))
 195         {
 196             indexOK = FALSE;
 197         }
 198         else
 199         {
 200             for (i = 0; i < ndon; i++)
 201             {
 202                 if (donindex[i] >= natoms)
 203                 {
 204                     indexOK = FALSE;
 205                 }
 206             }
 207             for (i = 0; i < nacc; i++)
 208             {
 209                 if (accindex[i] >= natoms)
 210                 {
 211                     indexOK = FALSE;
 212                 }
 213             }
 214         }
 215
 216         if (indexOK)
 217         {
 218
 219             if (bDatout)
 220             {
 221                 datfp = fopen(out_datfile, "w");
 222             }
 223
 224             if (bRKout)
 225             {
 226                 rkfp = xvgropen(out_xvgrkfile,
 227                                 "Distance and \\f{Symbol}k\\f{}\\S2\\N trajectory",
 228                                 "Time (ps)", "Distance (nm) / \\f{Symbol}k\\f{}\\S2\\N",
 229                                 oenv);
 230                 xvgr_legend(rkfp, 2, rkleg, oenv);
 231             }
 232
 233             if (bInstEffout)
 234             {
 235                 iefp = xvgropen(out_xvginstefffile,
 236                                 "Instantaneous RET Efficiency",
 237                                 "Time (ps)", "RET Efficiency",
 238                                 oenv);
 239                 xvgr_legend(iefp, 1, ieleg, oenv);
 240             }
 241
 242
 243             if (bRhistout)
 244             {
 245                 snew(rvalues, allocblock);
 246                 rblocksallocated += 1;
 247                 snew(rhist, histbins);
 248             }
 249
 250             if (bKhistout)
 251             {
 252                 snew(kappa2values, allocblock);
 253                 kblocksallocated += 1;
 254                 snew(khist, histbins);
 255             }
 256
 257             do
 258             {
 259                 clear_rvec(donvec);
 260                 clear_rvec(accvec);
 261                 clear_rvec(donpos);
 262                 clear_rvec(accpos);
 263                 for (i = 0; i < ndon / 2; i++)
 264                 {
 265                     rvec_sub(donvec, fr.x[donindex[2 * i]], donvec);
 266                     rvec_add(donvec, fr.x[donindex[2 * i + 1]], donvec);
 267                     rvec_add(donpos, fr.x[donindex[2 * i]], donpos);
 268                     rvec_add(donpos, fr.x[donindex[2 * i + 1]], donpos);
 269                 }
 270
 271                 for (i = 0; i < nacc / 2; i++)
 272                 {
 273                     rvec_sub(accvec, fr.x[accindex[2 * i]], accvec);
 274                     rvec_add(accvec, fr.x[accindex[2 * i + 1]], accvec);
 275                     rvec_add(accpos, fr.x[accindex[2 * i]], accpos);
 276                     rvec_add(accpos, fr.x[accindex[2 * i + 1]], accpos);
 277                 }
 278
 279                 unitv(donvec, donvec);
 280                 unitv(accvec, accvec);
 281
 282                 svmul((real) 1. / ndon, donpos, donpos);
 283                 svmul((real) 1. / nacc, accpos, accpos);
 284
 285                 if (bPBCdist)
 286                 {
 287                     set_pbc(pbc, ePBC, fr.box);
 288                     pbc_dx(pbc, donpos, accpos, dist);
 289                 }
 290                 else
 291                 {
 292                     rvec_sub(donpos, accpos, dist);
 293                 }
 294
 295                 unitv(dist, distnorm);
 296                 R       = norm(dist);
 297                 kappa2  = iprod(donvec, accvec)- 3.* (iprod(donvec, distnorm) * iprod(distnorm, accvec));
 298                 kappa2 *= kappa2;
 299                 if (R0 > 0)
 300                 {
 301                     Rfrac     = R/R0;
 302                     insteff   = 1/(1+(Rfrac*Rfrac*Rfrac*Rfrac*Rfrac*Rfrac)*2/3/kappa2);
 303                     insteffs += insteff;
 304
 305                     if (bInstEffout)
 306                     {
 307                         fprintf(iefp, "%12.7f %12.7f\n", fr.time, insteff);
 308                     }
 309                 }
 310
 311
 312                 Rs      += R;
 313                 kappa2s += kappa2;
 314                 rkcount++;
 315
 316                 if (bRKout)
 317                 {
 318                     fprintf(rkfp, "%12.7f %12.7f %12.7f\n", fr.time, R, kappa2);
 319                 }
 320
 321                 if (bDatout)
 322                 {
 323                     fprintf(datfp, "%12.7f %12.7f %12.7f\n", fr.time, R, kappa2);
 324                 }
 325
 326                 if (bRhistout)
 327                 {
 328                     rvalues[rkcount-1] = R;
 329                     if (rkcount % allocblock == 0)
 330                     {
 331                         srenew(rvalues, allocblock*(rblocksallocated+1));
 332                         rblocksallocated += 1;
 333                     }
 334                 }
 335
 336                 if (bKhistout)
 337                 {
 338                     kappa2values[rkcount-1] = kappa2;
 339                     if (rkcount % allocblock == 0)
 340                     {
 341                         srenew(kappa2values, allocblock*(kblocksallocated+1));
 342                         kblocksallocated += 1;
 343                     }
 344                 }
 345
 346                 bHaveNextFrame = read_next_frame(oenv, status, &fr);
 347             }
 348             while (bHaveNextFrame);
 349
 350             if (bRKout)
 351             {
 352                 ffclose(rkfp);
 353             }
 354
 355             if (bDatout)
 356             {
 357                 ffclose(datfp);
 358             }
 359
 360             if (bInstEffout)
 361             {
 362                 ffclose(iefp);
 363             }
 364
 365
 366             if (bRhistout)
 367             {
 368                 printf("Writing R-Histogram\n");
 369                 rmin = rvalues[0];
 370                 rmax = rvalues[0];
 371                 for (i = 1; i < rkcount; i++)
 372                 {
 373                     if (rvalues[i] < rmin)
 374                     {
 375                         rmin = rvalues[i];
 376                     }
 377                     else if (rvalues[i] > rmax)
 378                     {
 379                         rmax = rvalues[i];
 380                     }
 381                 }
 382                 rmin -= histexpand;
 383                 rmax += histexpand;
 384
 385                 rrange = rmax - rmin;
 386                 rincr  = rrange / histbins;
 387
 388                 for (i = 1; i < rkcount; i++)
 389                 {
 390                     bin         = (int) ((rvalues[i] - rmin) / rincr);
 391                     rhist[bin] += 1;
 392                 }
 393                 if (bNormHist)
 394                 {
 395                     for (i = 0; i < histbins; i++)
 396                     {
 397                         rhist[i] /= rkcount * rrange/histbins;
 398                     }
 399                     rhfp = xvgropen(out_xvgrhistfile, "Distance Distribution",
 400                                     "R (nm)", "Normalized Probability", oenv);
 401                 }
 402                 else
 403                 {
 404                     rhfp = xvgropen(out_xvgrhistfile, "Distance Distribution",
 405                                     "R (nm)", "Probability", oenv);
 406                 }
 407                 xvgr_legend(rhfp, 1, rhleg, oenv);
 408                 for (i = 0; i < histbins; i++)
 409                 {
 410                     fprintf(rhfp, "%12.7f %12.7f\n", (i + 0.5) * rincr + rmin,
 411                             rhist[i]);
 412                 }
 413                 ffclose(rhfp);
 414             }
 415
 416             if (bKhistout)
 417             {
 418                 printf("Writing kappa^2-Histogram\n");
 419                 krange = kmax - kmin;
 420                 kincr  = krange / histbins;
 421
 422                 for (i = 1; i < rkcount; i++)
 423                 {
 424                     bin         = (int) ((kappa2values[i] - kmin) / kincr);
 425                     khist[bin] += 1;
 426                 }
 427                 if (bNormHist)
 428                 {
 429                     for (i = 0; i < histbins; i++)
 430                     {
 431                         khist[i] /= rkcount * krange/histbins;
 432                     }
 433                     khfp = xvgropen(out_xvgkhistfile,
 434                                     "\\f{Symbol}k\\f{}\\S2\\N Distribution",
 435                                     "\\f{Symbol}k\\f{}\\S2\\N",
 436                                     "Normalized Probability", oenv);
 437                 }
 438                 else
 439                 {
 440                     khfp = xvgropen(out_xvgkhistfile,
 441                                     "\\f{Symbol}k\\f{}\\S2\\N Distribution",
 442                                     "\\f{Symbol}k\\f{}\\S2\\N", "Probability", oenv);
 443                 }
 444                 xvgr_legend(khfp, 1, khleg, oenv);
 445                 for (i = 0; i < histbins; i++)
 446                 {
 447                     fprintf(khfp, "%12.7f %12.7f\n", (i + 0.5) * kincr + kmin,
 448                             khist[i]);
 449                 }
 450                 ffclose(khfp);
 451             }
 452
 453             printf("\nAverages:\n");
 454             printf("R_avg   = %8.4f nm\nKappa^2 = %8.4f\n", Rs / rkcount,
 455                    kappa2s / rkcount);
 456             if (R0 > 0)
 457             {
 458                 printf("E_RETavg   = %8.4f\n", insteffs / rkcount);
 459             }
 460             please_cite(stdout, "Hoefling2011");
 461         }
 462         else
 463         {
 464             gmx_fatal(FARGS, "Index file invalid, check your index file for correct pairs.\n");
 465         }
 466     }
 467     else
 468     {
 469         gmx_fatal(FARGS, "Could not read first frame of the trajectory.\n");
 470     }
 471
 472     return 0;
 473 }