3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
9 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
10 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
11 * Copyright (c) 2001-2009, The GROMACS development team,
12 * check out http://www.gromacs.org for more information.
14 * This program is free software; you can redistribute it and/or
15 * modify it under the terms of the GNU General Public License
16 * as published by the Free Software Foundation; either version 2
17 * of the License, or (at your option) any later version.
19 * If you want to redistribute modifications, please consider that
20 * scientific software is very special. Version control is crucial -
21 * bugs must be traceable. We will be happy to consider code for
22 * inclusion in the official distribution, but derived work must not
23 * be called official GROMACS. Details are found in the README & COPYING
24 * files - if they are missing, get the official version at www.gromacs.org.
26 * To help us fund GROMACS development, we humbly ask that you cite
27 * the papers on the package - you can find them in the top README file.
29 * For more info, check our website at http://www.gromacs.org
41 int gmx_dyecoupl(int argc, char *argv[])
45 "This tool extracts dye dynamics from trajectory files.",
46 "Currently, R and kappa^2 between dyes is extracted for (F)RET",
47 "simulations with assumed dipolar coupling as in the Foerster equation.",
48 "It further allows the calculation of R(t) and kappa^2(t), R and",
49 "kappa^2 histograms and averages, as well as the instantaneous FRET",
50 "efficiency E(t) for a specified Foerster radius R_0 (switch [TT]-R0[tt]).",
51 "The input dyes have to be whole (see res and mol pbc options",
52 "in [TT]trjconv[tt]).",
53 "The dye transition dipole moment has to be defined by at least",
54 "a single atom pair, however multiple atom pairs can be provided ",
55 "in the index file. The distance R is calculated on the basis of",
56 "the COMs of the given atom pairs.",
57 "The [TT]-pbcdist[tt] option calculates distances to the nearest periodic",
58 "image instead to the distance in the box. This works however only,"
59 "for periodic boundaries in all 3 dimensions.",
60 "The [TT]-norm[tt] option (area-) normalizes the histograms."
63 static gmx_bool bPBCdist = FALSE, bNormHist = FALSE;
70 { "-pbcdist", FALSE, etBOOL, { &bPBCdist }, "Distance R based on PBC" },
71 { "-norm", FALSE, etBOOL, { &bNormHist }, "Normalize histograms" },
72 { "-bins", FALSE, etINT, {&histbins}, "# of histogram bins" },
73 { "-R0", FALSE, etREAL, {&R0}, "Foerster radius including kappa^2=2/3 in nm" }
79 { efTRX, "-f", NULL, ffREAD },
80 { efNDX, NULL, NULL, ffREAD },
81 { efXVG, "-ot", "rkappa", ffOPTWR },
82 { efXVG, "-oe", "insteff", ffOPTWR },
83 { efDAT, "-o", "rkappa", ffOPTWR },
84 { efXVG, "-rhist", "rhist", ffOPTWR },
85 { efXVG, "-khist", "khist", ffOPTWR }
87 #define NFILE asize(fnm)
90 const char *in_trajfile, *in_ndxfile, *out_xvgrkfile = NULL, *out_xvginstefffile = NULL, *out_xvgrhistfile = NULL, *out_xvgkhistfile = NULL, *out_datfile = NULL;
91 gmx_bool bHaveFirstFrame, bHaveNextFrame, indexOK = TRUE;
93 atom_id *donindex, *accindex;
95 t_atoms *atoms = NULL;
100 int allocblock = 1000;
101 real histexpand = 1e-6;
102 rvec donvec, accvec, donpos, accpos, dist, distnorm;
105 /*we rely on PBC autodetection (...currently)*/
108 real *rvalues = NULL, *kappa2values = NULL, *rhist = NULL, *khist = NULL;
111 FILE *rkfp = NULL, *rhfp = NULL, *khfp = NULL, *datfp = NULL, *iefp = NULL;
112 gmx_bool bRKout, bRhistout, bKhistout, bDatout, bInstEffout, grident;
114 const char *rkleg[2] = { "R", "\\f{Symbol}k\\f{}\\S2\\N" };
115 const char *rhleg[1] = { "p(R)" };
116 const char *khleg[1] = { "p(\\f{Symbol}k\\f{}\\S2\\N)" };
117 const char *ieleg[1] = { "E\\sRET\\N(t)" };
119 real R, kappa2, insteff, Rs = 0., kappa2s = 0., insteffs = 0., rmax, rmin, kmin = 0., kmax = 4.,
120 rrange, krange, rincr, kincr, Rfrac;
121 int rkcount = 0, rblocksallocated = 0, kblocksallocated = 0;
123 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);
126 /* Check command line options for filenames and set bool flags when switch used*/
127 in_trajfile = opt2fn("-f", NFILE, fnm);
128 in_ndxfile = opt2fn("-n", NFILE, fnm);
129 out_xvgrkfile = opt2fn("-ot", NFILE, fnm);
130 out_xvgrhistfile = opt2fn("-rhist", NFILE, fnm);
131 out_xvgkhistfile = opt2fn("-khist", NFILE, fnm);
132 out_xvginstefffile = opt2fn("-oe", NFILE, fnm);
133 out_datfile = opt2fn("-o", NFILE, fnm);
135 bRKout = opt2bSet("-ot", NFILE, fnm);
136 bRhistout = opt2bSet("-rhist", NFILE, fnm);
137 bKhistout = opt2bSet("-khist", NFILE, fnm);
138 bDatout = opt2bSet("-o", NFILE, fnm);
139 bInstEffout = opt2bSet("-oe", NFILE, fnm);
145 printf("Calculating distances to periodic image.\n");
146 printf("Be careful! This produces only valid results for PBC in all three dimensions\n");
150 if (bInstEffout && R0 <= 0.)
152 gmx_fatal(FARGS, "You have to specify R0 and R0 has to be larger than 0 nm.\n\n");
155 printf("Select group with donor atom pairs defining the transition moment\n");
156 get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &ndon, &donindex, &grpnm);
158 printf("Select group with acceptor atom pairs defining the transition moment\n");
159 get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &nacc, &accindex, &grpnm);
161 /*check if groups are identical*/
166 for (i = 0; i < nacc; i++)
168 if (accindex[i] != donindex[i])
178 gmx_fatal(FARGS, "Donor and acceptor group are identical. This makes no sense.");
181 printf("Reading first frame\n");
182 /* open trx file for reading */
184 flags = flags | TRX_READ_X;
185 bHaveFirstFrame = read_first_frame(oenv, &status, in_trajfile, &fr, flags);
189 printf("First frame is OK\n");
191 if ((ndon % 2 != 0) || (nacc % 2 != 0))
197 for (i = 0; i < ndon; i++)
199 if (donindex[i] >= natoms)
204 for (i = 0; i < nacc; i++)
206 if (accindex[i] >= natoms)
218 datfp = fopen(out_datfile, "w");
223 rkfp = xvgropen(out_xvgrkfile,
224 "Distance and \\f{Symbol}k\\f{}\\S2\\N trajectory",
225 "Time (ps)", "Distance (nm) / \\f{Symbol}k\\f{}\\S2\\N",
227 xvgr_legend(rkfp, 2, rkleg, oenv);
232 iefp = xvgropen(out_xvginstefffile,
233 "Instantaneous RET Efficiency",
234 "Time (ps)", "RET Efficiency",
236 xvgr_legend(iefp, 1, ieleg, oenv);
242 snew(rvalues, allocblock);
243 rblocksallocated += 1;
244 snew(rhist, histbins);
249 snew(kappa2values, allocblock);
250 kblocksallocated += 1;
251 snew(khist, histbins);
260 for (i = 0; i < ndon / 2; i++)
262 rvec_sub(donvec, fr.x[donindex[2 * i]], donvec);
263 rvec_add(donvec, fr.x[donindex[2 * i + 1]], donvec);
264 rvec_add(donpos, fr.x[donindex[2 * i]], donpos);
265 rvec_add(donpos, fr.x[donindex[2 * i + 1]], donpos);
268 for (i = 0; i < nacc / 2; i++)
270 rvec_sub(accvec, fr.x[accindex[2 * i]], accvec);
271 rvec_add(accvec, fr.x[accindex[2 * i + 1]], accvec);
272 rvec_add(accpos, fr.x[accindex[2 * i]], accpos);
273 rvec_add(accpos, fr.x[accindex[2 * i + 1]], accpos);
276 unitv(donvec, donvec);
277 unitv(accvec, accvec);
279 svmul((real) 1. / ndon, donpos, donpos);
280 svmul((real) 1. / nacc, accpos, accpos);
284 set_pbc(pbc, ePBC, fr.box);
285 pbc_dx(pbc, donpos, accpos, dist);
289 rvec_sub(donpos, accpos, dist);
292 unitv(dist, distnorm);
294 kappa2 = iprod(donvec, accvec)- 3.* (iprod(donvec, distnorm) * iprod(distnorm, accvec));
299 insteff = 1/(1+(Rfrac*Rfrac*Rfrac*Rfrac*Rfrac*Rfrac)*2/3/kappa2);
304 fprintf(iefp, "%12.7f %12.7f\n", fr.time, insteff);
315 fprintf(rkfp, "%12.7f %12.7f %12.7f\n", fr.time, R, kappa2);
320 fprintf(datfp, "%12.7f %12.7f %12.7f\n", fr.time, R, kappa2);
325 rvalues[rkcount-1] = R;
326 if (rkcount % allocblock == 0)
328 srenew(rvalues, allocblock*(rblocksallocated+1));
329 rblocksallocated += 1;
335 kappa2values[rkcount-1] = kappa2;
336 if (rkcount % allocblock == 0)
338 srenew(kappa2values, allocblock*(kblocksallocated+1));
339 kblocksallocated += 1;
343 bHaveNextFrame = read_next_frame(oenv, status, &fr);
345 while (bHaveNextFrame);
365 printf("Writing R-Histogram\n");
368 for (i = 1; i < rkcount; i++)
370 if (rvalues[i] < rmin)
374 else if (rvalues[i] > rmax)
382 rrange = rmax - rmin;
383 rincr = rrange / histbins;
385 for (i = 1; i < rkcount; i++)
387 bin = (int) ((rvalues[i] - rmin) / rincr);
392 for (i = 0; i < histbins; i++)
394 rhist[i] /= rkcount * rrange/histbins;
396 rhfp = xvgropen(out_xvgrhistfile, "Distance Distribution",
397 "R (nm)", "Normalized Probability", oenv);
401 rhfp = xvgropen(out_xvgrhistfile, "Distance Distribution",
402 "R (nm)", "Probability", oenv);
404 xvgr_legend(rhfp, 1, rhleg, oenv);
405 for (i = 0; i < histbins; i++)
407 fprintf(rhfp, "%12.7f %12.7f\n", (i + 0.5) * rincr + rmin,
415 printf("Writing kappa^2-Histogram\n");
416 krange = kmax - kmin;
417 kincr = krange / histbins;
419 for (i = 1; i < rkcount; i++)
421 bin = (int) ((kappa2values[i] - kmin) / kincr);
426 for (i = 0; i < histbins; i++)
428 khist[i] /= rkcount * krange/histbins;
430 khfp = xvgropen(out_xvgkhistfile,
431 "\\f{Symbol}k\\f{}\\S2\\N Distribution",
432 "\\f{Symbol}k\\f{}\\S2\\N",
433 "Normalized Probability", oenv);
437 khfp = xvgropen(out_xvgkhistfile,
438 "\\f{Symbol}k\\f{}\\S2\\N Distribution",
439 "\\f{Symbol}k\\f{}\\S2\\N", "Probability", oenv);
441 xvgr_legend(khfp, 1, khleg, oenv);
442 for (i = 0; i < histbins; i++)
444 fprintf(khfp, "%12.7f %12.7f\n", (i + 0.5) * kincr + kmin,
450 printf("\nAverages:\n");
451 printf("R_avg = %8.4f nm\nKappa^2 = %8.4f\n", Rs / rkcount,
455 printf("E_RETavg = %8.4f\n", insteffs / rkcount);
457 please_cite(stdout, "Hoefling2011");
461 gmx_fatal(FARGS, "Index file invalid, check your index file for correct pairs.\n");
466 gmx_fatal(FARGS, "Could not read first frame of the trajectory.\n");