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     PbcType pbcType = PbcType::Unset;
 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,
 129                            argv,
 130                            PCA_CAN_BEGIN | PCA_CAN_END | PCA_CAN_VIEW | PCA_TIME_UNIT,
 131                            NFILE,
 132                            fnm,
 133                            NPA,
 134                            pa,
 135                            asize(desc),
 136                            desc,
 137                            0,
 138                            nullptr,
 139                            &oenv))
 140     {
 141         return 0;
 142     }
 143
 144
 145     /* Check command line options for filenames and set bool flags when switch used*/
 146     in_trajfile        = opt2fn("-f", NFILE, fnm);
 147     out_xvgrkfile      = opt2fn("-ot", NFILE, fnm);
 148     out_xvgrhistfile   = opt2fn("-rhist", NFILE, fnm);
 149     out_xvgkhistfile   = opt2fn("-khist", NFILE, fnm);
 150     out_xvginstefffile = opt2fn("-oe", NFILE, fnm);
 151     out_datfile        = opt2fn("-o", NFILE, fnm);
 152
 153     bRKout      = opt2bSet("-ot", NFILE, fnm);
 154     bRhistout   = opt2bSet("-rhist", NFILE, fnm);
 155     bKhistout   = opt2bSet("-khist", NFILE, fnm);
 156     bDatout     = opt2bSet("-o", NFILE, fnm);
 157     bInstEffout = opt2bSet("-oe", NFILE, fnm);
 158
 159
 160     /* PBC warning. */
 161     if (bPBCdist)
 162     {
 163         printf("Calculating distances to periodic image.\n");
 164         printf("Be careful! This produces only valid results for PBC in all three dimensions\n");
 165     }
 166
 167
 168     if (bInstEffout && R0 <= 0.)
 169     {
 170         gmx_fatal(FARGS, "You have to specify R0 and R0 has to be larger than 0 nm.\n\n");
 171     }
 172
 173     printf("Select group with donor atom pairs defining the transition moment\n");
 174     get_index(nullptr, ftp2fn_null(efNDX, NFILE, fnm), 1, &ndon, &donindex, &grpnm);
 175
 176     printf("Select group with acceptor atom pairs defining the transition moment\n");
 177     get_index(nullptr, ftp2fn_null(efNDX, NFILE, fnm), 1, &nacc, &accindex, &grpnm);
 178
 179     /*check if groups are identical*/
 180     grident = TRUE;
 181
 182     if (ndon == nacc)
 183     {
 184         for (i = 0; i < nacc; i++)
 185         {
 186             if (accindex[i] != donindex[i])
 187             {
 188                 grident = FALSE;
 189                 break;
 190             }
 191         }
 192     }
 193
 194     if (grident)
 195     {
 196         gmx_fatal(FARGS, "Donor and acceptor group are identical. This makes no sense.");
 197     }
 198
 199     printf("Reading first frame\n");
 200     /* open trx file for reading */
 201     flags           = 0;
 202     flags           = flags | TRX_READ_X;
 203     bHaveFirstFrame = read_first_frame(oenv, &status, in_trajfile, &fr, flags);
 204
 205     if (bHaveFirstFrame)
 206     {
 207         printf("First frame is OK\n");
 208         natoms = fr.natoms;
 209         if ((ndon % 2 != 0) || (nacc % 2 != 0))
 210         {
 211             indexOK = FALSE;
 212         }
 213         else
 214         {
 215             for (i = 0; i < ndon; i++)
 216             {
 217                 if (donindex[i] >= natoms)
 218                 {
 219                     indexOK = FALSE;
 220                 }
 221             }
 222             for (i = 0; i < nacc; i++)
 223             {
 224                 if (accindex[i] >= natoms)
 225                 {
 226                     indexOK = FALSE;
 227                 }
 228             }
 229         }
 230
 231         if (indexOK)
 232         {
 233
 234             if (bDatout)
 235             {
 236                 datfp = gmx_ffopen(out_datfile, "w");
 237             }
 238
 239             if (bRKout)
 240             {
 241                 rkfp = xvgropen(out_xvgrkfile,
 242                                 "Distance and \\f{Symbol}k\\f{}\\S2\\N trajectory",
 243                                 "Time (ps)",
 244                                 "Distance (nm) / \\f{Symbol}k\\f{}\\S2\\N",
 245                                 oenv);
 246                 xvgr_legend(rkfp, 2, rkleg, oenv);
 247             }
 248
 249             if (bInstEffout)
 250             {
 251                 iefp = xvgropen(out_xvginstefffile,
 252                                 "Instantaneous RET Efficiency",
 253                                 "Time (ps)",
 254                                 "RET Efficiency",
 255                                 oenv);
 256                 xvgr_legend(iefp, 1, ieleg, oenv);
 257             }
 258
 259
 260             if (bRhistout)
 261             {
 262                 snew(rvalues, allocblock);
 263                 rblocksallocated += 1;
 264                 snew(rhist, histbins);
 265             }
 266
 267             if (bKhistout)
 268             {
 269                 snew(kappa2values, allocblock);
 270                 kblocksallocated += 1;
 271                 snew(khist, histbins);
 272             }
 273
 274             do
 275             {
 276                 clear_rvec(donvec);
 277                 clear_rvec(accvec);
 278                 clear_rvec(donpos);
 279                 clear_rvec(accpos);
 280                 for (i = 0; i < ndon / 2; i++)
 281                 {
 282                     rvec_sub(donvec, fr.x[donindex[2 * i]], donvec);
 283                     rvec_add(donvec, fr.x[donindex[2 * i + 1]], donvec);
 284                     rvec_add(donpos, fr.x[donindex[2 * i]], donpos);
 285                     rvec_add(donpos, fr.x[donindex[2 * i + 1]], donpos);
 286                 }
 287
 288                 for (i = 0; i < nacc / 2; i++)
 289                 {
 290                     rvec_sub(accvec, fr.x[accindex[2 * i]], accvec);
 291                     rvec_add(accvec, fr.x[accindex[2 * i + 1]], accvec);
 292                     rvec_add(accpos, fr.x[accindex[2 * i]], accpos);
 293                     rvec_add(accpos, fr.x[accindex[2 * i + 1]], accpos);
 294                 }
 295
 296                 unitv(donvec, donvec);
 297                 unitv(accvec, accvec);
 298
 299                 svmul(1.0 / ndon, donpos, donpos);
 300                 svmul(1.0 / nacc, accpos, accpos);
 301
 302                 if (bPBCdist)
 303                 {
 304                     set_pbc(pbc, pbcType, fr.box);
 305                     pbc_dx(pbc, donpos, accpos, dist);
 306                 }
 307                 else
 308                 {
 309                     rvec_sub(donpos, accpos, dist);
 310                 }
 311
 312                 unitv(dist, distnorm);
 313                 R = norm(dist);
 314                 kappa2 = iprod(donvec, accvec) - 3. * (iprod(donvec, distnorm) * iprod(distnorm, accvec));
 315                 kappa2 *= kappa2;
 316                 if (R0 > 0)
 317                 {
 318                     Rfrac = R / R0;
 319                     insteff = 1 / (1 + (Rfrac * Rfrac * Rfrac * Rfrac * Rfrac * Rfrac) * 2 / 3 / kappa2);
 320                     insteffs += insteff;
 321
 322                     if (bInstEffout)
 323                     {
 324                         fprintf(iefp, "%12.7f %12.7f\n", fr.time, insteff);
 325                     }
 326                 }
 327
 328
 329                 Rs += R;
 330                 kappa2s += kappa2;
 331                 rkcount++;
 332
 333                 if (bRKout)
 334                 {
 335                     fprintf(rkfp, "%12.7f %12.7f %12.7f\n", fr.time, R, kappa2);
 336                 }
 337
 338                 if (bDatout)
 339                 {
 340                     fprintf(datfp, "%12.7f %12.7f %12.7f\n", fr.time, R, kappa2);
 341                 }
 342
 343                 if (bRhistout)
 344                 {
 345                     rvalues[rkcount - 1] = R;
 346                     if (rkcount % allocblock == 0)
 347                     {
 348                         srenew(rvalues, allocblock * (rblocksallocated + 1));
 349                         rblocksallocated += 1;
 350                     }
 351                 }
 352
 353                 if (bKhistout)
 354                 {
 355                     kappa2values[rkcount - 1] = kappa2;
 356                     if (rkcount % allocblock == 0)
 357                     {
 358                         srenew(kappa2values, allocblock * (kblocksallocated + 1));
 359                         kblocksallocated += 1;
 360                     }
 361                 }
 362
 363                 bHaveNextFrame = read_next_frame(oenv, status, &fr);
 364             } while (bHaveNextFrame);
 365
 366             if (bRKout)
 367             {
 368                 xvgrclose(rkfp);
 369             }
 370
 371             if (bDatout)
 372             {
 373                 gmx_ffclose(datfp);
 374             }
 375
 376             if (bInstEffout)
 377             {
 378                 xvgrclose(iefp);
 379             }
 380
 381
 382             if (bRhistout)
 383             {
 384                 printf("Writing R-Histogram\n");
 385                 rmin = rvalues[0];
 386                 rmax = rvalues[0];
 387                 for (i = 1; i < rkcount; i++)
 388                 {
 389                     if (rvalues[i] < rmin)
 390                     {
 391                         rmin = rvalues[i];
 392                     }
 393                     else if (rvalues[i] > rmax)
 394                     {
 395                         rmax = rvalues[i];
 396                     }
 397                 }
 398                 rmin -= histexpand;
 399                 rmax += histexpand;
 400
 401                 rrange = rmax - rmin;
 402                 rincr  = rrange / histbins;
 403
 404                 for (i = 1; i < rkcount; i++)
 405                 {
 406                     bin = static_cast<int>((rvalues[i] - rmin) / rincr);
 407                     rhist[bin] += 1;
 408                 }
 409                 if (bNormHist)
 410                 {
 411                     for (i = 0; i < histbins; i++)
 412                     {
 413                         rhist[i] /= rkcount * rrange / histbins;
 414                     }
 415                     rhfp = xvgropen(out_xvgrhistfile,
 416                                     "Distance Distribution",
 417                                     "R (nm)",
 418                                     "Normalized Probability",
 419                                     oenv);
 420                 }
 421                 else
 422                 {
 423                     rhfp = xvgropen(
 424                             out_xvgrhistfile, "Distance Distribution", "R (nm)", "Probability", oenv);
 425                 }
 426                 xvgr_legend(rhfp, 1, rhleg, oenv);
 427                 for (i = 0; i < histbins; i++)
 428                 {
 429                     fprintf(rhfp, "%12.7f %12.7f\n", (i + 0.5) * rincr + rmin, rhist[i]);
 430                 }
 431                 xvgrclose(rhfp);
 432             }
 433
 434             if (bKhistout)
 435             {
 436                 printf("Writing kappa^2-Histogram\n");
 437                 krange = kmax - kmin;
 438                 kincr  = krange / histbins;
 439
 440                 for (i = 1; i < rkcount; i++)
 441                 {
 442                     bin = static_cast<int>((kappa2values[i] - kmin) / kincr);
 443                     khist[bin] += 1;
 444                 }
 445                 if (bNormHist)
 446                 {
 447                     for (i = 0; i < histbins; i++)
 448                     {
 449                         khist[i] /= rkcount * krange / histbins;
 450                     }
 451                     khfp = xvgropen(out_xvgkhistfile,
 452                                     "\\f{Symbol}k\\f{}\\S2\\N Distribution",
 453                                     "\\f{Symbol}k\\f{}\\S2\\N",
 454                                     "Normalized Probability",
 455                                     oenv);
 456                 }
 457                 else
 458                 {
 459                     khfp = xvgropen(out_xvgkhistfile,
 460                                     "\\f{Symbol}k\\f{}\\S2\\N Distribution",
 461                                     "\\f{Symbol}k\\f{}\\S2\\N",
 462                                     "Probability",
 463                                     oenv);
 464                 }
 465                 xvgr_legend(khfp, 1, khleg, oenv);
 466                 for (i = 0; i < histbins; i++)
 467                 {
 468                     fprintf(khfp, "%12.7f %12.7f\n", (i + 0.5) * kincr + kmin, khist[i]);
 469                 }
 470                 xvgrclose(khfp);
 471             }
 472
 473             printf("\nAverages:\n");
 474             printf("R_avg   = %8.4f nm\nKappa^2 = %8.4f\n", Rs / rkcount, kappa2s / rkcount);
 475             if (R0 > 0)
 476             {
 477                 printf("E_RETavg   = %8.4f\n", insteffs / rkcount);
 478             }
 479             please_cite(stdout, "Hoefling2011");
 480         }
 481         else
 482         {
 483             gmx_fatal(FARGS, "Index file invalid, check your index file for correct pairs.\n");
 484         }
 485     }
 486     else
 487     {
 488         gmx_fatal(FARGS, "Could not read first frame of the trajectory.\n");
 489     }
 490
 491     return 0;
 492 }