src/gromacs/gmxana/gmx_dyecoupl.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013,2014,2015,2017 by the GROMACS development team.
   5  * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   6  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   7  * and including many others, as listed in the AUTHORS file in the
   8  * top-level source directory and at http://www.gromacs.org.
   9  *
  10  * GROMACS is free software; you can redistribute it and/or
  11  * modify it under the terms of the GNU Lesser General Public License
  12  * as published by the Free Software Foundation; either version 2.1
  13  * of the License, or (at your option) any later version.
  14  *
  15  * GROMACS is distributed in the hope that it will be useful,
  16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  18  * Lesser General Public License for more details.
  19  *
  20  * You should have received a copy of the GNU Lesser General Public
  21  * License along with GROMACS; if not, see
  22  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  23  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  24  *
  25  * If you want to redistribute modifications to GROMACS, please
  26  * consider that scientific software is very special. Version
  27  * control is crucial - bugs must be traceable. We will be happy to
  28  * consider code for inclusion in the official distribution, but
  29  * derived work must not be called official GROMACS. Details are found
  30  * in the README & COPYING files - if they are missing, get the
  31  * official version at http://www.gromacs.org.
  32  *
  33  * To help us fund GROMACS development, we humbly ask that you cite
  34  * the research papers on the package. Check out http://www.gromacs.org.
  35  */
  36 #include "gmxpre.h"
  37
  38 #include "gromacs/commandline/filenm.h"
  39 #include "gromacs/commandline/pargs.h"
  40 #include "gromacs/fileio/trxio.h"
  41 #include "gromacs/fileio/xvgr.h"
  42 #include "gromacs/gmxana/gmx_ana.h"
  43 #include "gromacs/math/vec.h"
  44 #include "gromacs/pbcutil/pbc.h"
  45 #include "gromacs/topology/index.h"
  46 #include "gromacs/trajectory/trajectoryframe.h"
  47 #include "gromacs/utility/arraysize.h"
  48 #include "gromacs/utility/fatalerror.h"
  49 #include "gromacs/utility/futil.h"
  50 #include "gromacs/utility/pleasecite.h"
  51 #include "gromacs/utility/smalloc.h"
  52
  53 int gmx_dyecoupl(int argc, char* argv[])
  54 {
  55     const char* desc[] = {
  56         "[THISMODULE] extracts dye dynamics from trajectory files.",
  57         "Currently, R and kappa^2 between dyes is extracted for (F)RET",
  58         "simulations with assumed dipolar coupling as in the Foerster equation.",
  59         "It further allows the calculation of R(t) and kappa^2(t), R and",
  60         "kappa^2 histograms and averages, as well as the instantaneous FRET",
  61         "efficiency E(t) for a specified Foerster radius R_0 (switch [TT]-R0[tt]).",
  62         "The input dyes have to be whole (see res and mol pbc options",
  63         "in [TT]trjconv[tt]).",
  64         "The dye transition dipole moment has to be defined by at least",
  65         "a single atom pair, however multiple atom pairs can be provided ",
  66         "in the index file. The distance R is calculated on the basis of",
  67         "the COMs of the given atom pairs.",
  68         "The [TT]-pbcdist[tt] option calculates distances to the nearest periodic",
  69         "image instead to the distance in the box. This works however only,",
  70         "for periodic boundaries in all 3 dimensions.",
  71         "The [TT]-norm[tt] option (area-) normalizes the histograms."
  72     };
  73
  74     static gmx_bool   bPBCdist = FALSE, bNormHist = FALSE;
  75     int               histbins = 50;
  76     gmx_output_env_t* oenv;
  77     real              R0 = -1;
  78
  79     t_pargs pa[] = {
  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[] = { { efTRX, "-f", nullptr, ffREAD },     { efNDX, nullptr, nullptr, ffREAD },
  88                        { efXVG, "-ot", "rkappa", ffOPTWR },  { efXVG, "-oe", "insteff", ffOPTWR },
  89                        { efDAT, "-o", "rkappa", ffOPTWR },   { efXVG, "-rhist", "rhist", ffOPTWR },
  90                        { efXVG, "-khist", "khist", ffOPTWR } };
  91 #define NFILE asize(fnm)
  92
  93
  94     const char *in_trajfile, *out_xvgrkfile = nullptr, *out_xvginstefffile = nullptr,
  95                              *out_xvgrhistfile = nullptr, *out_xvgkhistfile = nullptr,
  96                              *out_datfile = nullptr;
  97     gmx_bool     bHaveFirstFrame, bHaveNextFrame, indexOK = TRUE;
  98     int          ndon, nacc;
  99     int *        donindex, *accindex;
 100     char*        grpnm;
 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 = nullptr, *kappa2values = nullptr, *rhist = nullptr, *khist = nullptr;
 114     t_pbc*   pbc = nullptr;
 115     int      i, bin;
 116     FILE *   rkfp = nullptr, *rhfp = nullptr, *khfp = nullptr, *datfp = nullptr, *iefp = nullptr;
 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,
 129                            NFILE, fnm, NPA, pa, asize(desc), desc, 0, nullptr, &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     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(nullptr, 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(nullptr, 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 = gmx_ffopen(out_datfile, "w");
 227             }
 228
 229             if (bRKout)
 230             {
 231                 rkfp = xvgropen(out_xvgrkfile, "Distance and \\f{Symbol}k\\f{}\\S2\\N trajectory",
 232                                 "Time (ps)", "Distance (nm) / \\f{Symbol}k\\f{}\\S2\\N", oenv);
 233                 xvgr_legend(rkfp, 2, rkleg, oenv);
 234             }
 235
 236             if (bInstEffout)
 237             {
 238                 iefp = xvgropen(out_xvginstefffile, "Instantaneous RET Efficiency", "Time (ps)",
 239                                 "RET Efficiency", 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(1.0 / ndon, donpos, donpos);
 284                 svmul(1.0 / 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             } while (bHaveNextFrame);
 349
 350             if (bRKout)
 351             {
 352                 xvgrclose(rkfp);
 353             }
 354
 355             if (bDatout)
 356             {
 357                 gmx_ffclose(datfp);
 358             }
 359
 360             if (bInstEffout)
 361             {
 362                 xvgrclose(iefp);
 363             }
 364
 365
 366             if (bRhistout)
 367             {
 368                 printf("Writing R-Histogram\n");
 369                 rmin = rvalues[0];
 370                 rmax = rvalues[0];
 371                 for (i = 1; i < rkcount; i++)
 372                 {
 373                     if (rvalues[i] < rmin)
 374                     {
 375                         rmin = rvalues[i];
 376                     }
 377                     else if (rvalues[i] > rmax)
 378                     {
 379                         rmax = rvalues[i];
 380                     }
 381                 }
 382                 rmin -= histexpand;
 383                 rmax += histexpand;
 384
 385                 rrange = rmax - rmin;
 386                 rincr  = rrange / histbins;
 387
 388                 for (i = 1; i < rkcount; i++)
 389                 {
 390                     bin = static_cast<int>((rvalues[i] - rmin) / rincr);
 391                     rhist[bin] += 1;
 392                 }
 393                 if (bNormHist)
 394                 {
 395                     for (i = 0; i < histbins; i++)
 396                     {
 397                         rhist[i] /= rkcount * rrange / histbins;
 398                     }
 399                     rhfp = xvgropen(out_xvgrhistfile, "Distance Distribution", "R (nm)",
 400                                     "Normalized Probability", oenv);
 401                 }
 402                 else
 403                 {
 404                     rhfp = xvgropen(out_xvgrhistfile, "Distance Distribution", "R (nm)",
 405                                     "Probability", oenv);
 406                 }
 407                 xvgr_legend(rhfp, 1, rhleg, oenv);
 408                 for (i = 0; i < histbins; i++)
 409                 {
 410                     fprintf(rhfp, "%12.7f %12.7f\n", (i + 0.5) * rincr + rmin, rhist[i]);
 411                 }
 412                 xvgrclose(rhfp);
 413             }
 414
 415             if (bKhistout)
 416             {
 417                 printf("Writing kappa^2-Histogram\n");
 418                 krange = kmax - kmin;
 419                 kincr  = krange / histbins;
 420
 421                 for (i = 1; i < rkcount; i++)
 422                 {
 423                     bin = static_cast<int>((kappa2values[i] - kmin) / kincr);
 424                     khist[bin] += 1;
 425                 }
 426                 if (bNormHist)
 427                 {
 428                     for (i = 0; i < histbins; i++)
 429                     {
 430                         khist[i] /= rkcount * krange / histbins;
 431                     }
 432                     khfp = xvgropen(out_xvgkhistfile, "\\f{Symbol}k\\f{}\\S2\\N Distribution",
 433                                     "\\f{Symbol}k\\f{}\\S2\\N", "Normalized Probability", oenv);
 434                 }
 435                 else
 436                 {
 437                     khfp = xvgropen(out_xvgkhistfile, "\\f{Symbol}k\\f{}\\S2\\N Distribution",
 438                                     "\\f{Symbol}k\\f{}\\S2\\N", "Probability", oenv);
 439                 }
 440                 xvgr_legend(khfp, 1, khleg, oenv);
 441                 for (i = 0; i < histbins; i++)
 442                 {
 443                     fprintf(khfp, "%12.7f %12.7f\n", (i + 0.5) * kincr + kmin, khist[i]);
 444                 }
 445                 xvgrclose(khfp);
 446             }
 447
 448             printf("\nAverages:\n");
 449             printf("R_avg   = %8.4f nm\nKappa^2 = %8.4f\n", Rs / rkcount, kappa2s / rkcount);
 450             if (R0 > 0)
 451             {
 452                 printf("E_RETavg   = %8.4f\n", insteffs / rkcount);
 453             }
 454             please_cite(stdout, "Hoefling2011");
 455         }
 456         else
 457         {
 458             gmx_fatal(FARGS, "Index file invalid, check your index file for correct pairs.\n");
 459         }
 460     }
 461     else
 462     {
 463         gmx_fatal(FARGS, "Could not read first frame of the trajectory.\n");
 464     }
 465
 466     return 0;
 467 }