2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
36 #include "gromacs/fileio/filenm.h"
38 #include "gromacs/math/vec.h"
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"
49 int gmx_dyecoupl(int argc, char *argv[])
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."
71 static gmx_bool bPBCdist = FALSE, bNormHist = FALSE;
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" }
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 }
95 #define NFILE asize(fnm)
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;
101 atom_id *donindex, *accindex;
107 int allocblock = 1000;
108 real histexpand = 1e-6;
109 rvec donvec, accvec, donpos, accpos, dist, distnorm;
112 /*we rely on PBC autodetection (...currently)*/
115 real *rvalues = NULL, *kappa2values = NULL, *rhist = NULL, *khist = NULL;
118 FILE *rkfp = NULL, *rhfp = NULL, *khfp = NULL, *datfp = NULL, *iefp = NULL;
119 gmx_bool bRKout, bRhistout, bKhistout, bDatout, bInstEffout, grident;
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)" };
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;
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))
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);
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);
155 printf("Calculating distances to periodic image.\n");
156 printf("Be careful! This produces only valid results for PBC in all three dimensions\n");
160 if (bInstEffout && R0 <= 0.)
162 gmx_fatal(FARGS, "You have to specify R0 and R0 has to be larger than 0 nm.\n\n");
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);
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);
171 /*check if groups are identical*/
176 for (i = 0; i < nacc; i++)
178 if (accindex[i] != donindex[i])
188 gmx_fatal(FARGS, "Donor and acceptor group are identical. This makes no sense.");
191 printf("Reading first frame\n");
192 /* open trx file for reading */
194 flags = flags | TRX_READ_X;
195 bHaveFirstFrame = read_first_frame(oenv, &status, in_trajfile, &fr, flags);
199 printf("First frame is OK\n");
201 if ((ndon % 2 != 0) || (nacc % 2 != 0))
207 for (i = 0; i < ndon; i++)
209 if (donindex[i] >= natoms)
214 for (i = 0; i < nacc; i++)
216 if (accindex[i] >= natoms)
228 datfp = fopen(out_datfile, "w");
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",
237 xvgr_legend(rkfp, 2, rkleg, oenv);
242 iefp = xvgropen(out_xvginstefffile,
243 "Instantaneous RET Efficiency",
244 "Time (ps)", "RET Efficiency",
246 xvgr_legend(iefp, 1, ieleg, oenv);
252 snew(rvalues, allocblock);
253 rblocksallocated += 1;
254 snew(rhist, histbins);
259 snew(kappa2values, allocblock);
260 kblocksallocated += 1;
261 snew(khist, histbins);
270 for (i = 0; i < ndon / 2; i++)
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);
278 for (i = 0; i < nacc / 2; i++)
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);
286 unitv(donvec, donvec);
287 unitv(accvec, accvec);
289 svmul((real) 1. / ndon, donpos, donpos);
290 svmul((real) 1. / nacc, accpos, accpos);
294 set_pbc(pbc, ePBC, fr.box);
295 pbc_dx(pbc, donpos, accpos, dist);
299 rvec_sub(donpos, accpos, dist);
302 unitv(dist, distnorm);
304 kappa2 = iprod(donvec, accvec)- 3.* (iprod(donvec, distnorm) * iprod(distnorm, accvec));
309 insteff = 1/(1+(Rfrac*Rfrac*Rfrac*Rfrac*Rfrac*Rfrac)*2/3/kappa2);
314 fprintf(iefp, "%12.7f %12.7f\n", fr.time, insteff);
325 fprintf(rkfp, "%12.7f %12.7f %12.7f\n", fr.time, R, kappa2);
330 fprintf(datfp, "%12.7f %12.7f %12.7f\n", fr.time, R, kappa2);
335 rvalues[rkcount-1] = R;
336 if (rkcount % allocblock == 0)
338 srenew(rvalues, allocblock*(rblocksallocated+1));
339 rblocksallocated += 1;
345 kappa2values[rkcount-1] = kappa2;
346 if (rkcount % allocblock == 0)
348 srenew(kappa2values, allocblock*(kblocksallocated+1));
349 kblocksallocated += 1;
353 bHaveNextFrame = read_next_frame(oenv, status, &fr);
355 while (bHaveNextFrame);
375 printf("Writing R-Histogram\n");
378 for (i = 1; i < rkcount; i++)
380 if (rvalues[i] < rmin)
384 else if (rvalues[i] > rmax)
392 rrange = rmax - rmin;
393 rincr = rrange / histbins;
395 for (i = 1; i < rkcount; i++)
397 bin = (int) ((rvalues[i] - rmin) / rincr);
402 for (i = 0; i < histbins; i++)
404 rhist[i] /= rkcount * rrange/histbins;
406 rhfp = xvgropen(out_xvgrhistfile, "Distance Distribution",
407 "R (nm)", "Normalized Probability", oenv);
411 rhfp = xvgropen(out_xvgrhistfile, "Distance Distribution",
412 "R (nm)", "Probability", oenv);
414 xvgr_legend(rhfp, 1, rhleg, oenv);
415 for (i = 0; i < histbins; i++)
417 fprintf(rhfp, "%12.7f %12.7f\n", (i + 0.5) * rincr + rmin,
425 printf("Writing kappa^2-Histogram\n");
426 krange = kmax - kmin;
427 kincr = krange / histbins;
429 for (i = 1; i < rkcount; i++)
431 bin = (int) ((kappa2values[i] - kmin) / kincr);
436 for (i = 0; i < histbins; i++)
438 khist[i] /= rkcount * krange/histbins;
440 khfp = xvgropen(out_xvgkhistfile,
441 "\\f{Symbol}k\\f{}\\S2\\N Distribution",
442 "\\f{Symbol}k\\f{}\\S2\\N",
443 "Normalized Probability", oenv);
447 khfp = xvgropen(out_xvgkhistfile,
448 "\\f{Symbol}k\\f{}\\S2\\N Distribution",
449 "\\f{Symbol}k\\f{}\\S2\\N", "Probability", oenv);
451 xvgr_legend(khfp, 1, khleg, oenv);
452 for (i = 0; i < histbins; i++)
454 fprintf(khfp, "%12.7f %12.7f\n", (i + 0.5) * kincr + kmin,
460 printf("\nAverages:\n");
461 printf("R_avg = %8.4f nm\nKappa^2 = %8.4f\n", Rs / rkcount,
465 printf("E_RETavg = %8.4f\n", insteffs / rkcount);
467 please_cite(stdout, "Hoefling2011");
471 gmx_fatal(FARGS, "Index file invalid, check your index file for correct pairs.\n");
476 gmx_fatal(FARGS, "Could not read first frame of the trajectory.\n");