2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2009, The GROMACS development team,
6 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012,2013, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
52 int gmx_dyecoupl(int argc, char *argv[])
56 "This tool extracts dye dynamics from trajectory files.",
57 "Currently, R and kappa^2 between dyes is extracted for (F)RET",
58 "simulations with assumed dipolar coupling as in the Foerster equation.",
59 "It further allows the calculation of R(t) and kappa^2(t), R and",
60 "kappa^2 histograms and averages, as well as the instantaneous FRET",
61 "efficiency E(t) for a specified Foerster radius R_0 (switch [TT]-R0[tt]).",
62 "The input dyes have to be whole (see res and mol pbc options",
63 "in [TT]trjconv[tt]).",
64 "The dye transition dipole moment has to be defined by at least",
65 "a single atom pair, however multiple atom pairs can be provided ",
66 "in the index file. The distance R is calculated on the basis of",
67 "the COMs of the given atom pairs.",
68 "The [TT]-pbcdist[tt] option calculates distances to the nearest periodic",
69 "image instead to the distance in the box. This works however only,"
70 "for periodic boundaries in all 3 dimensions.",
71 "The [TT]-norm[tt] option (area-) normalizes the histograms."
74 static gmx_bool bPBCdist = FALSE, bNormHist = FALSE;
81 { "-pbcdist", FALSE, etBOOL, { &bPBCdist },"Distance R based on PBC" },
82 { "-norm", FALSE, etBOOL, { &bNormHist },"Normalize histograms" },
83 { "-bins", FALSE, etINT, {&histbins},"# of histogram bins" },
84 { "-R0", FALSE, etREAL, {&R0},"Foerster radius including kappa^2=2/3 in nm" }
90 { efTRX, "-f", NULL, ffREAD },
91 { efNDX, NULL, NULL, ffREAD },
92 { efXVG, "-ot", "rkappa",ffOPTWR },
93 { efXVG, "-oe", "insteff",ffOPTWR },
94 { efDAT, "-o", "rkappa",ffOPTWR },
95 { efXVG, "-rhist","rhist", ffOPTWR },
96 { efXVG, "-khist", "khist", ffOPTWR }
98 #define NFILE asize(fnm)
101 const char *in_trajfile, *in_ndxfile, *out_xvgrkfile = NULL, *out_xvginstefffile = NULL, *out_xvgrhistfile = NULL, *out_xvgkhistfile = NULL,*out_datfile=NULL;
102 gmx_bool bHaveFirstFrame, bHaveNextFrame, indexOK = TRUE;
104 atom_id *donindex, *accindex;
106 t_atoms *atoms = NULL;
111 int allocblock = 1000;
112 real histexpand = 1e-6;
113 rvec donvec, accvec, donpos, accpos, dist, distnorm;
116 /*we rely on PBC autodetection (...currently)*/
119 real *rvalues=NULL, *kappa2values=NULL, *rhist=NULL, *khist=NULL;
122 FILE *rkfp = NULL, *rhfp = NULL, *khfp = NULL,*datfp=NULL,*iefp=NULL;
123 gmx_bool bRKout, bRhistout, bKhistout,bDatout,bInstEffout,grident;
125 const char *rkleg[2] = { "R", "\\f{Symbol}k\\f{}\\S2\\N" };
126 const char *rhleg[1] = { "p(R)" };
127 const char *khleg[1] = { "p(\\f{Symbol}k\\f{}\\S2\\N)" };
128 const char *ieleg[1] = { "E\\sRET\\N(t)" };
130 real R, kappa2, insteff, Rs = 0., kappa2s = 0., insteffs=0., rmax, rmin, kmin = 0., kmax = 4.,
131 rrange, krange, rincr, kincr,Rfrac;
132 int rkcount = 0, rblocksallocated = 0, kblocksallocated = 0;
134 CopyRight(stderr, argv[0]);
135 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);
138 /* Check command line options for filenames and set bool flags when switch used*/
139 in_trajfile = opt2fn("-f", NFILE, fnm);
140 in_ndxfile = opt2fn("-n", NFILE, fnm);
141 out_xvgrkfile = opt2fn("-ot", NFILE, fnm);
142 out_xvgrhistfile = opt2fn("-rhist", NFILE, fnm);
143 out_xvgkhistfile = opt2fn("-khist", NFILE, fnm);
144 out_xvginstefffile = opt2fn("-oe", NFILE, fnm);
145 out_datfile = opt2fn("-o",NFILE,fnm);
147 bRKout = opt2bSet("-ot", NFILE, fnm);
148 bRhistout = opt2bSet("-rhist", NFILE, fnm);
149 bKhistout = opt2bSet("-khist", NFILE, fnm);
150 bDatout = opt2bSet("-o", NFILE, fnm);
151 bInstEffout = opt2bSet("-oe", NFILE, fnm);
157 printf("Calculating distances to periodic image.\n");
158 printf("Be careful! This produces only valid results for PBC in all three dimensions\n");
162 if (bInstEffout && R0<=0.)
164 gmx_fatal(FARGS,"You have to specify R0 and R0 has to be larger than 0 nm.\n\n");
167 printf("Select group with donor atom pairs defining the transition moment\n");
168 get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &ndon, &donindex,&grpnm);
170 printf("Select group with acceptor atom pairs defining the transition moment\n");
171 get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &nacc, &accindex,&grpnm);
173 /*check if groups are identical*/
180 if(accindex[i] != donindex[i])
190 gmx_fatal(FARGS,"Donor and acceptor group are identical. This makes no sense.");
193 printf("Reading first frame\n");
194 /* open trx file for reading */
196 flags = flags | TRX_READ_X;
197 bHaveFirstFrame = read_first_frame(oenv, &status, in_trajfile, &fr, flags);
201 printf("First frame is OK\n");
203 if ((ndon % 2 != 0) || (nacc % 2 != 0))
209 for (i = 0; i < ndon;i++)
211 if (donindex[i] >= natoms)
214 for (i = 0; i < nacc;i++)
216 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);
322 fprintf(rkfp, "%12.7f %12.7f %12.7f\n", fr.time, R, kappa2);
325 fprintf(datfp, "%12.7f %12.7f %12.7f\n", fr.time, R, kappa2);
329 rvalues[rkcount-1] = R;
330 if (rkcount % allocblock == 0)
332 srenew(rvalues, allocblock*(rblocksallocated+1));
333 rblocksallocated += 1;
339 kappa2values[rkcount-1] = kappa2;
340 if (rkcount % allocblock == 0)
342 srenew(kappa2values, allocblock*(kblocksallocated+1));
343 kblocksallocated += 1;
347 bHaveNextFrame = read_next_frame(oenv, status, &fr);
348 } while (bHaveNextFrame);
362 printf("Writing R-Histogram\n");
365 for (i = 1; i < rkcount; i++)
367 if (rvalues[i] < rmin)
369 else if (rvalues[i] > rmax)
375 rrange = rmax - rmin;
376 rincr = rrange / histbins;
378 for (i = 1; i < rkcount; i++)
380 bin = (int) ((rvalues[i] - rmin) / rincr);
385 for (i = 0; i < histbins; i++)
386 rhist[i] /= rkcount * rrange/histbins;
387 rhfp = xvgropen(out_xvgrhistfile, "Distance Distribution",
388 "R (nm)", "Normalized Probability", oenv);
391 rhfp = xvgropen(out_xvgrhistfile, "Distance Distribution",
392 "R (nm)", "Probability", oenv);
394 xvgr_legend(rhfp, 1, rhleg, oenv);
395 for (i = 0; i < histbins; i++)
397 fprintf(rhfp, "%12.7f %12.7f\n", (i + 0.5) * rincr + rmin,
405 printf("Writing kappa^2-Histogram\n");
406 krange = kmax - kmin;
407 kincr = krange / histbins;
409 for (i = 1; i < rkcount; i++)
411 bin = (int) ((kappa2values[i] - kmin) / kincr);
416 for (i = 0; i < histbins; i++)
417 khist[i] /= rkcount * krange/histbins;
418 khfp = xvgropen(out_xvgkhistfile,
419 "\\f{Symbol}k\\f{}\\S2\\N Distribution",
420 "\\f{Symbol}k\\f{}\\S2\\N",
421 "Normalized Probability", oenv);
424 khfp = xvgropen(out_xvgkhistfile,
425 "\\f{Symbol}k\\f{}\\S2\\N Distribution",
426 "\\f{Symbol}k\\f{}\\S2\\N", "Probability", oenv);
428 xvgr_legend(khfp, 1, khleg, oenv);
429 for (i = 0; i < histbins; i++)
431 fprintf(khfp, "%12.7f %12.7f\n", (i + 0.5) * kincr + kmin,
437 printf("\nAverages:\n");
438 printf("R_avg = %8.4f nm\nKappa^2 = %8.4f\n", Rs / rkcount,
442 printf("E_RETavg = %8.4f\n", insteffs / rkcount);
444 please_cite(stdout,"Hoefling2011");
448 gmx_fatal(FARGS,"Index file invalid, check your index file for correct pairs.\n");
453 gmx_fatal(FARGS,"Could not read first frame of the trajectory.\n");