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