src/gromacs/gmxana/gmx_dyecoupl.c

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