2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
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"
53 int gmx_dyecoupl(int argc, char* argv[])
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."
74 static gmx_bool bPBCdist = FALSE, bNormHist = FALSE;
76 gmx_output_env_t* oenv;
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" }
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)
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;
99 int * donindex, *accindex;
105 int allocblock = 1000;
106 real histexpand = 1e-6;
107 rvec donvec, accvec, donpos, accpos, dist, distnorm;
110 /*we rely on PBC autodetection (...currently)*/
111 PbcType pbcType = PbcType::Unset;
113 real * rvalues = nullptr, *kappa2values = nullptr, *rhist = nullptr, *khist = nullptr;
114 t_pbc* pbc = nullptr;
116 FILE * rkfp = nullptr, *rhfp = nullptr, *khfp = nullptr, *datfp = nullptr, *iefp = nullptr;
117 gmx_bool bRKout, bRhistout, bKhistout, bDatout, bInstEffout, grident;
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)" };
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;
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))
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);
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);
153 printf("Calculating distances to periodic image.\n");
154 printf("Be careful! This produces only valid results for PBC in all three dimensions\n");
158 if (bInstEffout && R0 <= 0.)
160 gmx_fatal(FARGS, "You have to specify R0 and R0 has to be larger than 0 nm.\n\n");
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);
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);
169 /*check if groups are identical*/
174 for (i = 0; i < nacc; i++)
176 if (accindex[i] != donindex[i])
186 gmx_fatal(FARGS, "Donor and acceptor group are identical. This makes no sense.");
189 printf("Reading first frame\n");
190 /* open trx file for reading */
192 flags = flags | TRX_READ_X;
193 bHaveFirstFrame = read_first_frame(oenv, &status, in_trajfile, &fr, flags);
197 printf("First frame is OK\n");
199 if ((ndon % 2 != 0) || (nacc % 2 != 0))
205 for (i = 0; i < ndon; i++)
207 if (donindex[i] >= natoms)
212 for (i = 0; i < nacc; i++)
214 if (accindex[i] >= natoms)
226 datfp = gmx_ffopen(out_datfile, "w");
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);
238 iefp = xvgropen(out_xvginstefffile, "Instantaneous RET Efficiency", "Time (ps)",
239 "RET Efficiency", oenv);
240 xvgr_legend(iefp, 1, ieleg, oenv);
246 snew(rvalues, allocblock);
247 rblocksallocated += 1;
248 snew(rhist, histbins);
253 snew(kappa2values, allocblock);
254 kblocksallocated += 1;
255 snew(khist, histbins);
264 for (i = 0; i < ndon / 2; i++)
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);
272 for (i = 0; i < nacc / 2; i++)
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);
280 unitv(donvec, donvec);
281 unitv(accvec, accvec);
283 svmul(1.0 / ndon, donpos, donpos);
284 svmul(1.0 / nacc, accpos, accpos);
288 set_pbc(pbc, pbcType, fr.box);
289 pbc_dx(pbc, donpos, accpos, dist);
293 rvec_sub(donpos, accpos, dist);
296 unitv(dist, distnorm);
298 kappa2 = iprod(donvec, accvec) - 3. * (iprod(donvec, distnorm) * iprod(distnorm, accvec));
303 insteff = 1 / (1 + (Rfrac * Rfrac * Rfrac * Rfrac * Rfrac * Rfrac) * 2 / 3 / kappa2);
308 fprintf(iefp, "%12.7f %12.7f\n", fr.time, insteff);
319 fprintf(rkfp, "%12.7f %12.7f %12.7f\n", fr.time, R, kappa2);
324 fprintf(datfp, "%12.7f %12.7f %12.7f\n", fr.time, R, kappa2);
329 rvalues[rkcount - 1] = R;
330 if (rkcount % allocblock == 0)
332 srenew(rvalues, allocblock * (rblocksallocated + 1));
333 rblocksallocated += 1;
339 kappa2values[rkcount - 1] = kappa2;
340 if (rkcount % allocblock == 0)
342 srenew(kappa2values, allocblock * (kblocksallocated + 1));
343 kblocksallocated += 1;
347 bHaveNextFrame = read_next_frame(oenv, status, &fr);
348 } while (bHaveNextFrame);
368 printf("Writing R-Histogram\n");
371 for (i = 1; i < rkcount; i++)
373 if (rvalues[i] < rmin)
377 else if (rvalues[i] > rmax)
385 rrange = rmax - rmin;
386 rincr = rrange / histbins;
388 for (i = 1; i < rkcount; i++)
390 bin = static_cast<int>((rvalues[i] - rmin) / rincr);
395 for (i = 0; i < histbins; i++)
397 rhist[i] /= rkcount * rrange / histbins;
399 rhfp = xvgropen(out_xvgrhistfile, "Distance Distribution", "R (nm)",
400 "Normalized Probability", oenv);
404 rhfp = xvgropen(out_xvgrhistfile, "Distance Distribution", "R (nm)",
405 "Probability", oenv);
407 xvgr_legend(rhfp, 1, rhleg, oenv);
408 for (i = 0; i < histbins; i++)
410 fprintf(rhfp, "%12.7f %12.7f\n", (i + 0.5) * rincr + rmin, rhist[i]);
417 printf("Writing kappa^2-Histogram\n");
418 krange = kmax - kmin;
419 kincr = krange / histbins;
421 for (i = 1; i < rkcount; i++)
423 bin = static_cast<int>((kappa2values[i] - kmin) / kincr);
428 for (i = 0; i < histbins; i++)
430 khist[i] /= rkcount * krange / histbins;
432 khfp = xvgropen(out_xvgkhistfile, "\\f{Symbol}k\\f{}\\S2\\N Distribution",
433 "\\f{Symbol}k\\f{}\\S2\\N", "Normalized Probability", oenv);
437 khfp = xvgropen(out_xvgkhistfile, "\\f{Symbol}k\\f{}\\S2\\N Distribution",
438 "\\f{Symbol}k\\f{}\\S2\\N", "Probability", oenv);
440 xvgr_legend(khfp, 1, khleg, oenv);
441 for (i = 0; i < histbins; i++)
443 fprintf(khfp, "%12.7f %12.7f\n", (i + 0.5) * kincr + kmin, khist[i]);
448 printf("\nAverages:\n");
449 printf("R_avg = %8.4f nm\nKappa^2 = %8.4f\n", Rs / rkcount, kappa2s / rkcount);
452 printf("E_RETavg = %8.4f\n", insteffs / rkcount);
454 please_cite(stdout, "Hoefling2011");
458 gmx_fatal(FARGS, "Index file invalid, check your index file for correct pairs.\n");
463 gmx_fatal(FARGS, "Could not read first frame of the trajectory.\n");