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