2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2012,2013, 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"
40 #include "gromacs/commandline/pargs.h"
43 #include "gromacs/fileio/trxio.h"
46 int gmx_dyecoupl(int argc, char *argv[])
50 "[THISMODULE] extracts dye dynamics from trajectory files.",
51 "Currently, R and kappa^2 between dyes is extracted for (F)RET",
52 "simulations with assumed dipolar coupling as in the Foerster equation.",
53 "It further allows the calculation of R(t) and kappa^2(t), R and",
54 "kappa^2 histograms and averages, as well as the instantaneous FRET",
55 "efficiency E(t) for a specified Foerster radius R_0 (switch [TT]-R0[tt]).",
56 "The input dyes have to be whole (see res and mol pbc options",
57 "in [TT]trjconv[tt]).",
58 "The dye transition dipole moment has to be defined by at least",
59 "a single atom pair, however multiple atom pairs can be provided ",
60 "in the index file. The distance R is calculated on the basis of",
61 "the COMs of the given atom pairs.",
62 "The [TT]-pbcdist[tt] option calculates distances to the nearest periodic",
63 "image instead to the distance in the box. This works however only,"
64 "for periodic boundaries in all 3 dimensions.",
65 "The [TT]-norm[tt] option (area-) normalizes the histograms."
68 static gmx_bool bPBCdist = FALSE, bNormHist = FALSE;
75 { "-pbcdist", FALSE, etBOOL, { &bPBCdist }, "Distance R based on PBC" },
76 { "-norm", FALSE, etBOOL, { &bNormHist }, "Normalize histograms" },
77 { "-bins", FALSE, etINT, {&histbins}, "# of histogram bins" },
78 { "-R0", FALSE, etREAL, {&R0}, "Foerster radius including kappa^2=2/3 in nm" }
84 { efTRX, "-f", NULL, ffREAD },
85 { efNDX, NULL, NULL, ffREAD },
86 { efXVG, "-ot", "rkappa", ffOPTWR },
87 { efXVG, "-oe", "insteff", ffOPTWR },
88 { efDAT, "-o", "rkappa", ffOPTWR },
89 { efXVG, "-rhist", "rhist", ffOPTWR },
90 { efXVG, "-khist", "khist", ffOPTWR }
92 #define NFILE asize(fnm)
95 const char *in_trajfile, *in_ndxfile, *out_xvgrkfile = NULL, *out_xvginstefffile = NULL, *out_xvgrhistfile = NULL, *out_xvgkhistfile = NULL, *out_datfile = NULL;
96 gmx_bool bHaveFirstFrame, bHaveNextFrame, indexOK = TRUE;
98 atom_id *donindex, *accindex;
100 t_atoms *atoms = NULL;
105 int allocblock = 1000;
106 real histexpand = 1e-6;
107 rvec donvec, accvec, donpos, accpos, dist, distnorm;
110 /*we rely on PBC autodetection (...currently)*/
113 real *rvalues = NULL, *kappa2values = NULL, *rhist = NULL, *khist = NULL;
116 FILE *rkfp = NULL, *rhfp = NULL, *khfp = NULL, *datfp = NULL, *iefp = NULL;
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 | PCA_BE_NICE, NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv))
134 /* Check command line options for filenames and set bool flags when switch used*/
135 in_trajfile = opt2fn("-f", NFILE, fnm);
136 in_ndxfile = opt2fn("-n", 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(atoms, 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(atoms, 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 = fopen(out_datfile, "w");
231 rkfp = xvgropen(out_xvgrkfile,
232 "Distance and \\f{Symbol}k\\f{}\\S2\\N trajectory",
233 "Time (ps)", "Distance (nm) / \\f{Symbol}k\\f{}\\S2\\N",
235 xvgr_legend(rkfp, 2, rkleg, oenv);
240 iefp = xvgropen(out_xvginstefffile,
241 "Instantaneous RET Efficiency",
242 "Time (ps)", "RET Efficiency",
244 xvgr_legend(iefp, 1, ieleg, oenv);
250 snew(rvalues, allocblock);
251 rblocksallocated += 1;
252 snew(rhist, histbins);
257 snew(kappa2values, allocblock);
258 kblocksallocated += 1;
259 snew(khist, histbins);
268 for (i = 0; i < ndon / 2; i++)
270 rvec_sub(donvec, fr.x[donindex[2 * i]], donvec);
271 rvec_add(donvec, fr.x[donindex[2 * i + 1]], donvec);
272 rvec_add(donpos, fr.x[donindex[2 * i]], donpos);
273 rvec_add(donpos, fr.x[donindex[2 * i + 1]], donpos);
276 for (i = 0; i < nacc / 2; i++)
278 rvec_sub(accvec, fr.x[accindex[2 * i]], accvec);
279 rvec_add(accvec, fr.x[accindex[2 * i + 1]], accvec);
280 rvec_add(accpos, fr.x[accindex[2 * i]], accpos);
281 rvec_add(accpos, fr.x[accindex[2 * i + 1]], accpos);
284 unitv(donvec, donvec);
285 unitv(accvec, accvec);
287 svmul((real) 1. / ndon, donpos, donpos);
288 svmul((real) 1. / nacc, accpos, accpos);
292 set_pbc(pbc, ePBC, fr.box);
293 pbc_dx(pbc, donpos, accpos, dist);
297 rvec_sub(donpos, accpos, dist);
300 unitv(dist, distnorm);
302 kappa2 = iprod(donvec, accvec)- 3.* (iprod(donvec, distnorm) * iprod(distnorm, accvec));
307 insteff = 1/(1+(Rfrac*Rfrac*Rfrac*Rfrac*Rfrac*Rfrac)*2/3/kappa2);
312 fprintf(iefp, "%12.7f %12.7f\n", fr.time, insteff);
323 fprintf(rkfp, "%12.7f %12.7f %12.7f\n", fr.time, R, kappa2);
328 fprintf(datfp, "%12.7f %12.7f %12.7f\n", fr.time, R, kappa2);
333 rvalues[rkcount-1] = R;
334 if (rkcount % allocblock == 0)
336 srenew(rvalues, allocblock*(rblocksallocated+1));
337 rblocksallocated += 1;
343 kappa2values[rkcount-1] = kappa2;
344 if (rkcount % allocblock == 0)
346 srenew(kappa2values, allocblock*(kblocksallocated+1));
347 kblocksallocated += 1;
351 bHaveNextFrame = read_next_frame(oenv, status, &fr);
353 while (bHaveNextFrame);
373 printf("Writing R-Histogram\n");
376 for (i = 1; i < rkcount; i++)
378 if (rvalues[i] < rmin)
382 else if (rvalues[i] > rmax)
390 rrange = rmax - rmin;
391 rincr = rrange / histbins;
393 for (i = 1; i < rkcount; i++)
395 bin = (int) ((rvalues[i] - rmin) / rincr);
400 for (i = 0; i < histbins; i++)
402 rhist[i] /= rkcount * rrange/histbins;
404 rhfp = xvgropen(out_xvgrhistfile, "Distance Distribution",
405 "R (nm)", "Normalized Probability", oenv);
409 rhfp = xvgropen(out_xvgrhistfile, "Distance Distribution",
410 "R (nm)", "Probability", oenv);
412 xvgr_legend(rhfp, 1, rhleg, oenv);
413 for (i = 0; i < histbins; i++)
415 fprintf(rhfp, "%12.7f %12.7f\n", (i + 0.5) * rincr + rmin,
423 printf("Writing kappa^2-Histogram\n");
424 krange = kmax - kmin;
425 kincr = krange / histbins;
427 for (i = 1; i < rkcount; i++)
429 bin = (int) ((kappa2values[i] - kmin) / kincr);
434 for (i = 0; i < histbins; i++)
436 khist[i] /= rkcount * krange/histbins;
438 khfp = xvgropen(out_xvgkhistfile,
439 "\\f{Symbol}k\\f{}\\S2\\N Distribution",
440 "\\f{Symbol}k\\f{}\\S2\\N",
441 "Normalized Probability", oenv);
445 khfp = xvgropen(out_xvgkhistfile,
446 "\\f{Symbol}k\\f{}\\S2\\N Distribution",
447 "\\f{Symbol}k\\f{}\\S2\\N", "Probability", oenv);
449 xvgr_legend(khfp, 1, khleg, oenv);
450 for (i = 0; i < histbins; i++)
452 fprintf(khfp, "%12.7f %12.7f\n", (i + 0.5) * kincr + kmin,
458 printf("\nAverages:\n");
459 printf("R_avg = %8.4f nm\nKappa^2 = %8.4f\n", Rs / rkcount,
463 printf("E_RETavg = %8.4f\n", insteffs / rkcount);
465 please_cite(stdout, "Hoefling2011");
469 gmx_fatal(FARGS, "Index file invalid, check your index file for correct pairs.\n");
474 gmx_fatal(FARGS, "Could not read first frame of the trajectory.\n");