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