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