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 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))
129 /* Check command line options for filenames and set bool flags when switch used*/
130 in_trajfile = opt2fn("-f", NFILE, fnm);
131 in_ndxfile = opt2fn("-n", NFILE, fnm);
132 out_xvgrkfile = opt2fn("-ot", NFILE, fnm);
133 out_xvgrhistfile = opt2fn("-rhist", NFILE, fnm);
134 out_xvgkhistfile = opt2fn("-khist", NFILE, fnm);
135 out_xvginstefffile = opt2fn("-oe", NFILE, fnm);
136 out_datfile = opt2fn("-o", NFILE, fnm);
138 bRKout = opt2bSet("-ot", NFILE, fnm);
139 bRhistout = opt2bSet("-rhist", NFILE, fnm);
140 bKhistout = opt2bSet("-khist", NFILE, fnm);
141 bDatout = opt2bSet("-o", NFILE, fnm);
142 bInstEffout = opt2bSet("-oe", NFILE, fnm);
148 printf("Calculating distances to periodic image.\n");
149 printf("Be careful! This produces only valid results for PBC in all three dimensions\n");
153 if (bInstEffout && R0 <= 0.)
155 gmx_fatal(FARGS, "You have to specify R0 and R0 has to be larger than 0 nm.\n\n");
158 printf("Select group with donor atom pairs defining the transition moment\n");
159 get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &ndon, &donindex, &grpnm);
161 printf("Select group with acceptor atom pairs defining the transition moment\n");
162 get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &nacc, &accindex, &grpnm);
164 /*check if groups are identical*/
169 for (i = 0; i < nacc; i++)
171 if (accindex[i] != donindex[i])
181 gmx_fatal(FARGS, "Donor and acceptor group are identical. This makes no sense.");
184 printf("Reading first frame\n");
185 /* open trx file for reading */
187 flags = flags | TRX_READ_X;
188 bHaveFirstFrame = read_first_frame(oenv, &status, in_trajfile, &fr, flags);
192 printf("First frame is OK\n");
194 if ((ndon % 2 != 0) || (nacc % 2 != 0))
200 for (i = 0; i < ndon; i++)
202 if (donindex[i] >= natoms)
207 for (i = 0; i < nacc; i++)
209 if (accindex[i] >= natoms)
221 datfp = fopen(out_datfile, "w");
226 rkfp = xvgropen(out_xvgrkfile,
227 "Distance and \\f{Symbol}k\\f{}\\S2\\N trajectory",
228 "Time (ps)", "Distance (nm) / \\f{Symbol}k\\f{}\\S2\\N",
230 xvgr_legend(rkfp, 2, rkleg, oenv);
235 iefp = xvgropen(out_xvginstefffile,
236 "Instantaneous RET Efficiency",
237 "Time (ps)", "RET Efficiency",
239 xvgr_legend(iefp, 1, ieleg, oenv);
245 snew(rvalues, allocblock);
246 rblocksallocated += 1;
247 snew(rhist, histbins);
252 snew(kappa2values, allocblock);
253 kblocksallocated += 1;
254 snew(khist, histbins);
263 for (i = 0; i < ndon / 2; i++)
265 rvec_sub(donvec, fr.x[donindex[2 * i]], donvec);
266 rvec_add(donvec, fr.x[donindex[2 * i + 1]], donvec);
267 rvec_add(donpos, fr.x[donindex[2 * i]], donpos);
268 rvec_add(donpos, fr.x[donindex[2 * i + 1]], donpos);
271 for (i = 0; i < nacc / 2; i++)
273 rvec_sub(accvec, fr.x[accindex[2 * i]], accvec);
274 rvec_add(accvec, fr.x[accindex[2 * i + 1]], accvec);
275 rvec_add(accpos, fr.x[accindex[2 * i]], accpos);
276 rvec_add(accpos, fr.x[accindex[2 * i + 1]], accpos);
279 unitv(donvec, donvec);
280 unitv(accvec, accvec);
282 svmul((real) 1. / ndon, donpos, donpos);
283 svmul((real) 1. / nacc, accpos, accpos);
287 set_pbc(pbc, ePBC, fr.box);
288 pbc_dx(pbc, donpos, accpos, dist);
292 rvec_sub(donpos, accpos, dist);
295 unitv(dist, distnorm);
297 kappa2 = iprod(donvec, accvec)- 3.* (iprod(donvec, distnorm) * iprod(distnorm, accvec));
302 insteff = 1/(1+(Rfrac*Rfrac*Rfrac*Rfrac*Rfrac*Rfrac)*2/3/kappa2);
307 fprintf(iefp, "%12.7f %12.7f\n", fr.time, insteff);
318 fprintf(rkfp, "%12.7f %12.7f %12.7f\n", fr.time, R, kappa2);
323 fprintf(datfp, "%12.7f %12.7f %12.7f\n", fr.time, R, kappa2);
328 rvalues[rkcount-1] = R;
329 if (rkcount % allocblock == 0)
331 srenew(rvalues, allocblock*(rblocksallocated+1));
332 rblocksallocated += 1;
338 kappa2values[rkcount-1] = kappa2;
339 if (rkcount % allocblock == 0)
341 srenew(kappa2values, allocblock*(kblocksallocated+1));
342 kblocksallocated += 1;
346 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 = (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",
400 "R (nm)", "Normalized Probability", oenv);
404 rhfp = xvgropen(out_xvgrhistfile, "Distance Distribution",
405 "R (nm)", "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,
418 printf("Writing kappa^2-Histogram\n");
419 krange = kmax - kmin;
420 kincr = krange / histbins;
422 for (i = 1; i < rkcount; i++)
424 bin = (int) ((kappa2values[i] - kmin) / kincr);
429 for (i = 0; i < histbins; i++)
431 khist[i] /= rkcount * krange/histbins;
433 khfp = xvgropen(out_xvgkhistfile,
434 "\\f{Symbol}k\\f{}\\S2\\N Distribution",
435 "\\f{Symbol}k\\f{}\\S2\\N",
436 "Normalized Probability", oenv);
440 khfp = xvgropen(out_xvgkhistfile,
441 "\\f{Symbol}k\\f{}\\S2\\N Distribution",
442 "\\f{Symbol}k\\f{}\\S2\\N", "Probability", oenv);
444 xvgr_legend(khfp, 1, khleg, oenv);
445 for (i = 0; i < histbins; i++)
447 fprintf(khfp, "%12.7f %12.7f\n", (i + 0.5) * kincr + kmin,
453 printf("\nAverages:\n");
454 printf("R_avg = %8.4f nm\nKappa^2 = %8.4f\n", Rs / rkcount,
458 printf("E_RETavg = %8.4f\n", insteffs / rkcount);
460 please_cite(stdout, "Hoefling2011");
464 gmx_fatal(FARGS, "Index file invalid, check your index file for correct pairs.\n");
469 gmx_fatal(FARGS, "Could not read first frame of the trajectory.\n");