src/gromacs/gmxana/gmx_dyecoupl.c

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