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