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