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"
39 #include "gromacs/utility/smalloc.h"
40 #include "gromacs/commandline/pargs.h"
41 #include "gromacs/math/vec.h"
42 #include "gromacs/fileio/xvgr.h"
43 #include "gromacs/fileio/trxio.h"
45 #include "gromacs/utility/fatalerror.h"
47 int gmx_dyecoupl(int argc, char *argv[])
51 "[THISMODULE] extracts dye dynamics from trajectory files.",
52 "Currently, R and kappa^2 between dyes is extracted for (F)RET",
53 "simulations with assumed dipolar coupling as in the Foerster equation.",
54 "It further allows the calculation of R(t) and kappa^2(t), R and",
55 "kappa^2 histograms and averages, as well as the instantaneous FRET",
56 "efficiency E(t) for a specified Foerster radius R_0 (switch [TT]-R0[tt]).",
57 "The input dyes have to be whole (see res and mol pbc options",
58 "in [TT]trjconv[tt]).",
59 "The dye transition dipole moment has to be defined by at least",
60 "a single atom pair, however multiple atom pairs can be provided ",
61 "in the index file. The distance R is calculated on the basis of",
62 "the COMs of the given atom pairs.",
63 "The [TT]-pbcdist[tt] option calculates distances to the nearest periodic",
64 "image instead to the distance in the box. This works however only,"
65 "for periodic boundaries in all 3 dimensions.",
66 "The [TT]-norm[tt] option (area-) normalizes the histograms."
69 static gmx_bool bPBCdist = FALSE, bNormHist = FALSE;
76 { "-pbcdist", FALSE, etBOOL, { &bPBCdist }, "Distance R based on PBC" },
77 { "-norm", FALSE, etBOOL, { &bNormHist }, "Normalize histograms" },
78 { "-bins", FALSE, etINT, {&histbins}, "# of histogram bins" },
79 { "-R0", FALSE, etREAL, {&R0}, "Foerster radius including kappa^2=2/3 in nm" }
85 { efTRX, "-f", NULL, ffREAD },
86 { efNDX, NULL, NULL, ffREAD },
87 { efXVG, "-ot", "rkappa", ffOPTWR },
88 { efXVG, "-oe", "insteff", ffOPTWR },
89 { efDAT, "-o", "rkappa", ffOPTWR },
90 { efXVG, "-rhist", "rhist", ffOPTWR },
91 { efXVG, "-khist", "khist", ffOPTWR }
93 #define NFILE asize(fnm)
96 const char *in_trajfile, *in_ndxfile, *out_xvgrkfile = NULL, *out_xvginstefffile = NULL, *out_xvgrhistfile = NULL, *out_xvgkhistfile = NULL, *out_datfile = NULL;
97 gmx_bool bHaveFirstFrame, bHaveNextFrame, indexOK = TRUE;
99 atom_id *donindex, *accindex;
101 t_atoms *atoms = NULL;
106 int allocblock = 1000;
107 real histexpand = 1e-6;
108 rvec donvec, accvec, donpos, accpos, dist, distnorm;
111 /*we rely on PBC autodetection (...currently)*/
114 real *rvalues = NULL, *kappa2values = NULL, *rhist = NULL, *khist = NULL;
117 FILE *rkfp = NULL, *rhfp = NULL, *khfp = NULL, *datfp = NULL, *iefp = NULL;
118 gmx_bool bRKout, bRhistout, bKhistout, bDatout, bInstEffout, grident;
120 const char *rkleg[2] = { "R", "\\f{Symbol}k\\f{}\\S2\\N" };
121 const char *rhleg[1] = { "p(R)" };
122 const char *khleg[1] = { "p(\\f{Symbol}k\\f{}\\S2\\N)" };
123 const char *ieleg[1] = { "E\\sRET\\N(t)" };
125 real R, kappa2, insteff, Rs = 0., kappa2s = 0., insteffs = 0., rmax, rmin, kmin = 0., kmax = 4.,
126 rrange, krange, rincr, kincr, Rfrac;
127 int rkcount = 0, rblocksallocated = 0, kblocksallocated = 0;
129 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))
135 /* Check command line options for filenames and set bool flags when switch used*/
136 in_trajfile = opt2fn("-f", NFILE, fnm);
137 in_ndxfile = opt2fn("-n", NFILE, fnm);
138 out_xvgrkfile = opt2fn("-ot", NFILE, fnm);
139 out_xvgrhistfile = opt2fn("-rhist", NFILE, fnm);
140 out_xvgkhistfile = opt2fn("-khist", NFILE, fnm);
141 out_xvginstefffile = opt2fn("-oe", NFILE, fnm);
142 out_datfile = opt2fn("-o", NFILE, fnm);
144 bRKout = opt2bSet("-ot", NFILE, fnm);
145 bRhistout = opt2bSet("-rhist", NFILE, fnm);
146 bKhistout = opt2bSet("-khist", NFILE, fnm);
147 bDatout = opt2bSet("-o", NFILE, fnm);
148 bInstEffout = opt2bSet("-oe", NFILE, fnm);
154 printf("Calculating distances to periodic image.\n");
155 printf("Be careful! This produces only valid results for PBC in all three dimensions\n");
159 if (bInstEffout && R0 <= 0.)
161 gmx_fatal(FARGS, "You have to specify R0 and R0 has to be larger than 0 nm.\n\n");
164 printf("Select group with donor atom pairs defining the transition moment\n");
165 get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &ndon, &donindex, &grpnm);
167 printf("Select group with acceptor atom pairs defining the transition moment\n");
168 get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &nacc, &accindex, &grpnm);
170 /*check if groups are identical*/
175 for (i = 0; i < nacc; i++)
177 if (accindex[i] != donindex[i])
187 gmx_fatal(FARGS, "Donor and acceptor group are identical. This makes no sense.");
190 printf("Reading first frame\n");
191 /* open trx file for reading */
193 flags = flags | TRX_READ_X;
194 bHaveFirstFrame = read_first_frame(oenv, &status, in_trajfile, &fr, flags);
198 printf("First frame is OK\n");
200 if ((ndon % 2 != 0) || (nacc % 2 != 0))
206 for (i = 0; i < ndon; i++)
208 if (donindex[i] >= natoms)
213 for (i = 0; i < nacc; i++)
215 if (accindex[i] >= natoms)
227 datfp = fopen(out_datfile, "w");
232 rkfp = xvgropen(out_xvgrkfile,
233 "Distance and \\f{Symbol}k\\f{}\\S2\\N trajectory",
234 "Time (ps)", "Distance (nm) / \\f{Symbol}k\\f{}\\S2\\N",
236 xvgr_legend(rkfp, 2, rkleg, oenv);
241 iefp = xvgropen(out_xvginstefffile,
242 "Instantaneous RET Efficiency",
243 "Time (ps)", "RET Efficiency",
245 xvgr_legend(iefp, 1, ieleg, oenv);
251 snew(rvalues, allocblock);
252 rblocksallocated += 1;
253 snew(rhist, histbins);
258 snew(kappa2values, allocblock);
259 kblocksallocated += 1;
260 snew(khist, histbins);
269 for (i = 0; i < ndon / 2; i++)
271 rvec_sub(donvec, fr.x[donindex[2 * i]], donvec);
272 rvec_add(donvec, fr.x[donindex[2 * i + 1]], donvec);
273 rvec_add(donpos, fr.x[donindex[2 * i]], donpos);
274 rvec_add(donpos, fr.x[donindex[2 * i + 1]], donpos);
277 for (i = 0; i < nacc / 2; i++)
279 rvec_sub(accvec, fr.x[accindex[2 * i]], accvec);
280 rvec_add(accvec, fr.x[accindex[2 * i + 1]], accvec);
281 rvec_add(accpos, fr.x[accindex[2 * i]], accpos);
282 rvec_add(accpos, fr.x[accindex[2 * i + 1]], accpos);
285 unitv(donvec, donvec);
286 unitv(accvec, accvec);
288 svmul((real) 1. / ndon, donpos, donpos);
289 svmul((real) 1. / nacc, accpos, accpos);
293 set_pbc(pbc, ePBC, fr.box);
294 pbc_dx(pbc, donpos, accpos, dist);
298 rvec_sub(donpos, accpos, dist);
301 unitv(dist, distnorm);
303 kappa2 = iprod(donvec, accvec)- 3.* (iprod(donvec, distnorm) * iprod(distnorm, accvec));
308 insteff = 1/(1+(Rfrac*Rfrac*Rfrac*Rfrac*Rfrac*Rfrac)*2/3/kappa2);
313 fprintf(iefp, "%12.7f %12.7f\n", fr.time, insteff);
324 fprintf(rkfp, "%12.7f %12.7f %12.7f\n", fr.time, R, kappa2);
329 fprintf(datfp, "%12.7f %12.7f %12.7f\n", fr.time, R, kappa2);
334 rvalues[rkcount-1] = R;
335 if (rkcount % allocblock == 0)
337 srenew(rvalues, allocblock*(rblocksallocated+1));
338 rblocksallocated += 1;
344 kappa2values[rkcount-1] = kappa2;
345 if (rkcount % allocblock == 0)
347 srenew(kappa2values, allocblock*(kblocksallocated+1));
348 kblocksallocated += 1;
352 bHaveNextFrame = read_next_frame(oenv, status, &fr);
354 while (bHaveNextFrame);
374 printf("Writing R-Histogram\n");
377 for (i = 1; i < rkcount; i++)
379 if (rvalues[i] < rmin)
383 else if (rvalues[i] > rmax)
391 rrange = rmax - rmin;
392 rincr = rrange / histbins;
394 for (i = 1; i < rkcount; i++)
396 bin = (int) ((rvalues[i] - rmin) / rincr);
401 for (i = 0; i < histbins; i++)
403 rhist[i] /= rkcount * rrange/histbins;
405 rhfp = xvgropen(out_xvgrhistfile, "Distance Distribution",
406 "R (nm)", "Normalized Probability", oenv);
410 rhfp = xvgropen(out_xvgrhistfile, "Distance Distribution",
411 "R (nm)", "Probability", oenv);
413 xvgr_legend(rhfp, 1, rhleg, oenv);
414 for (i = 0; i < histbins; i++)
416 fprintf(rhfp, "%12.7f %12.7f\n", (i + 0.5) * rincr + rmin,
424 printf("Writing kappa^2-Histogram\n");
425 krange = kmax - kmin;
426 kincr = krange / histbins;
428 for (i = 1; i < rkcount; i++)
430 bin = (int) ((kappa2values[i] - kmin) / kincr);
435 for (i = 0; i < histbins; i++)
437 khist[i] /= rkcount * krange/histbins;
439 khfp = xvgropen(out_xvgkhistfile,
440 "\\f{Symbol}k\\f{}\\S2\\N Distribution",
441 "\\f{Symbol}k\\f{}\\S2\\N",
442 "Normalized Probability", oenv);
446 khfp = xvgropen(out_xvgkhistfile,
447 "\\f{Symbol}k\\f{}\\S2\\N Distribution",
448 "\\f{Symbol}k\\f{}\\S2\\N", "Probability", oenv);
450 xvgr_legend(khfp, 1, khleg, oenv);
451 for (i = 0; i < histbins; i++)
453 fprintf(khfp, "%12.7f %12.7f\n", (i + 0.5) * kincr + kmin,
459 printf("\nAverages:\n");
460 printf("R_avg = %8.4f nm\nKappa^2 = %8.4f\n", Rs / rkcount,
464 printf("E_RETavg = %8.4f\n", insteffs / rkcount);
466 please_cite(stdout, "Hoefling2011");
470 gmx_fatal(FARGS, "Index file invalid, check your index file for correct pairs.\n");
475 gmx_fatal(FARGS, "Could not read first frame of the trajectory.\n");