From b3681c00b7ba591296b8c1fcfdc115e979c530d6 Mon Sep 17 00:00:00 2001 From: Carsten Kutzner Date: Tue, 20 Dec 2011 18:53:32 +0100 Subject: [PATCH] New tool that extracts dye dynamics from trajectories Change-Id: Ie2007a60ad612938160197c0ec091e474e7091d2 --- include/gmx_ana.h | 3 + src/gmxlib/copyrite.c | 8 +- src/tools/CMakeLists.txt | 4 +- src/tools/g_dyecoupl.c | 50 +++++ src/tools/gmx_dyecoupl.c | 431 +++++++++++++++++++++++++++++++++++++++ 5 files changed, 493 insertions(+), 3 deletions(-) create mode 100644 src/tools/g_dyecoupl.c create mode 100644 src/tools/gmx_dyecoupl.c diff --git a/include/gmx_ana.h b/include/gmx_ana.h index ecd659bd8f..cb885f068b 100644 --- a/include/gmx_ana.h +++ b/include/gmx_ana.h @@ -99,6 +99,9 @@ gmx_dist(int argc,char *argv[]); int gmx_dos(int argc,char *argv[]); +int +gmx_dyecoupl(int argc,char *argv[]); + int gmx_dyndom(int argc,char *argv[]); diff --git a/src/gmxlib/copyrite.c b/src/gmxlib/copyrite.c index 9b0d0285a3..606157f5b9 100644 --- a/src/gmxlib/copyrite.c +++ b/src/gmxlib/copyrite.c @@ -542,7 +542,13 @@ void please_cite(FILE *fp,const char *key) "C. Kutzner and J. Czub and H. Grubmuller", "Keep it Flexible: Driving Macromolecular Rotary Motions in Atomistic Simulations with GROMACS", "J. Chem. Theory Comput.", - 7, 2011, "1381-1393" } + 7, 2011, "1381-1393" }, + { "Hoefling2011", + "M. Hoefling, N. Lima, D. Haenni, C.A.M. Seidel, B. Schuler, H. Grubmuller", + "Structural Heterogeneity and Quantitative FRET Efficiency Distributions of Polyprolines through a Hybrid Atomistic Simulation and Monte Carlo Approach", + "PLoS ONE", + 6, 2011, "e19791" + } }; #define NSTR (int)asize(citedb) diff --git a/src/tools/CMakeLists.txt b/src/tools/CMakeLists.txt index 81d4f63449..ac3894d8b2 100644 --- a/src/tools/CMakeLists.txt +++ b/src/tools/CMakeLists.txt @@ -31,7 +31,7 @@ add_library(gmxana addconf.c calcpot.c edittop.c gmx_bar.c gmx_membed.c gmx_pme_error.c gmx_options.c gmx_dos.c gmx_hydorder.c gmx_densorder.c powerspect.c dens_filter.c - binsearch.c + binsearch.c gmx_dyecoupl.c ) @@ -55,7 +55,7 @@ set(GMX_TOOLS_PROGRAMS g_spol g_spatial g_tcaf g_traj g_tune_pme g_vanhove g_velacc g_clustsize g_mdmat g_wham g_sigeps g_bar g_membed g_pme_error g_rmsdist g_rotmat g_options - g_dos g_hydorder g_densorder + g_dos g_hydorder g_densorder g_dyecoupl ) set(GMX_TOOLS_PROGRAMS_NOT_FOR_INSTALLATION diff --git a/src/tools/g_dyecoupl.c b/src/tools/g_dyecoupl.c new file mode 100644 index 0000000000..73a5d329e0 --- /dev/null +++ b/src/tools/g_dyecoupl.c @@ -0,0 +1,50 @@ +/* + * + * This source code is part of + * + * G R O M A C S + * + * GROningen MAchine for Chemical Simulations + * + * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others. + * Copyright (c) 1991-2000, University of Groningen, The Netherlands. + * Copyright (c) 2001-2008, The GROMACS development team, + * check out http://www.gromacs.org for more information. + + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. + * + * If you want to redistribute modifications, please consider that + * scientific software is very special. Version control is crucial - + * bugs must be traceable. We will be happy to consider code for + * inclusion in the official distribution, but derived work must not + * be called official GROMACS. Details are found in the README & COPYING + * files - if they are missing, get the official version at www.gromacs.org. + * + * To help us fund GROMACS development, we humbly ask that you cite + * the papers on the package - you can find them in the top README file. + * + * For more info, check our website at http://www.gromacs.org + * + * And Hey: + * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include + + +/* This is just a wrapper binary. +*/ +int +main(int argc, char *argv[]) +{ + gmx_dyecoupl(argc,argv); + return 0; +} + diff --git a/src/tools/gmx_dyecoupl.c b/src/tools/gmx_dyecoupl.c new file mode 100644 index 0000000000..1817a5e5fd --- /dev/null +++ b/src/tools/gmx_dyecoupl.c @@ -0,0 +1,431 @@ +/* + * + * This source code is part of + * + * G R O M A C S + * + * GROningen MAchine for Chemical Simulations + * + * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others. + * Copyright (c) 1991-2000, University of Groningen, The Netherlands. + * Copyright (c) 2001-2009, The GROMACS development team, + * check out http://www.gromacs.org for more information. + + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. + * + * If you want to redistribute modifications, please consider that + * scientific software is very special. Version control is crucial - + * bugs must be traceable. We will be happy to consider code for + * inclusion in the official distribution, but derived work must not + * be called official GROMACS. Details are found in the README & COPYING + * files - if they are missing, get the official version at www.gromacs.org. + * + * To help us fund GROMACS development, we humbly ask that you cite + * the papers on the package - you can find them in the top README file. + * + * For more info, check our website at http://www.gromacs.org + */ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +int gmx_dyecoupl(int argc, char *argv[]) +{ + const char *desc[] = + { + "This tool extracts dye dynamics from trajectory files.", + "Currently, R and kappa^2 between dyes is extracted for (F)RET", + "simulations with assumed dipolar coupling as in the Foerster equation.", + "It further allows the calculation of R(t) and kappa^2(t), R and", + "kappa^2 histograms and averages, as well as the instantaneous FRET", + "efficiency E(t) for a specified Foerster radius R_0 (switch [TT]-R0[tt]).", + "The input dyes have to be whole (see res and mol pbc options", + "in [TT]trjconv[tt]).", + "The dye transition dipole moment has to be defined by at least", + "a single atom pair, however multiple atom pairs can be provided ", + "in the index file. The distance R is calculated on the basis of", + "the COMs of the given atom pairs.", + "The [TT]-pbcdist[tt] option calculates distances to the nearest periodic", + "image instead to the distance in the box. This works however only," + "for periodic boundaries in all 3 dimensions.", + "The [TT]-norm[tt] option (area-) normalizes the histograms." + }; + + static gmx_bool bPBCdist = FALSE, bNormHist = FALSE; + int histbins = 50; + output_env_t oenv; + real R0=-1; + + t_pargs pa[] = + { + { "-pbcdist", FALSE, etBOOL, { &bPBCdist },"Distance R based on PBC" }, + { "-norm", FALSE, etBOOL, { &bNormHist },"Normalize histograms" }, + { "-bins", FALSE, etINT, {&histbins},"# of histogram bins" }, + { "-R0", FALSE, etREAL, {&R0},"Foerster radius including kappa^2=2/3 in nm" } + }; +#define NPA asize(pa) + + t_filenm fnm[] = + { + { efTRX, "-f", NULL, ffREAD }, + { efNDX, NULL, NULL, ffREAD }, + { efXVG, "-ot", "rkappa",ffOPTWR }, + { efXVG, "-oe", "insteff",ffOPTWR }, + { efDAT, "-o", "rkappa",ffOPTWR }, + { efXVG, "-rhist","rhist", ffOPTWR }, + { efXVG, "-khist", "khist", ffOPTWR } + }; +#define NFILE asize(fnm) + + + const char *in_trajfile, *in_ndxfile, *out_xvgrkfile = NULL, *out_xvginstefffile = NULL, *out_xvgrhistfile = NULL, *out_xvgkhistfile = NULL,*out_datfile=NULL; + gmx_bool bHaveFirstFrame, bHaveNextFrame, indexOK = TRUE; + int ndon, nacc; + atom_id *donindex, *accindex; + char *grpnm; + t_atoms *atoms = NULL; + t_trxstatus *status; + t_trxframe fr; + + int flags; + int allocblock = 1000; + real histexpand = 1e-6; + rvec donvec, accvec, donpos, accpos, dist, distnorm; + int natoms; + + /*we rely on PBC autodetection (...currently)*/ + int ePBC = -1; + + real *rvalues=NULL, *kappa2values=NULL, *rhist=NULL, *khist=NULL; + t_pbc *pbc=NULL; + int i, bin; + FILE *rkfp = NULL, *rhfp = NULL, *khfp = NULL,*datfp=NULL,*iefp=NULL; + gmx_bool bRKout, bRhistout, bKhistout,bDatout,bInstEffout; + + const char *rkleg[2] = { "R", "\\f{Symbol}k\\f{}\\S2\\N" }; + const char *rhleg[1] = { "p(R)" }; + const char *khleg[1] = { "p(\\f{Symbol}k\\f{}\\S2\\N)" }; + const char *ieleg[1] = { "E\\sRET\\N(t)" }; + + real R, kappa2, insteff, Rs = 0., kappa2s = 0., insteffs=0., rmax, rmin, kmin = 0., kmax = 4., + rrange, krange, rincr, kincr,Rfrac; + int rkcount = 0, rblocksallocated = 0, kblocksallocated = 0; + + CopyRight(stderr, argv[0]); + 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); + + + /* Check command line options for filenames and set bool flags when switch used*/ + in_trajfile = opt2fn("-f", NFILE, fnm); + in_ndxfile = opt2fn("-n", NFILE, fnm); + out_xvgrkfile = opt2fn("-ot", NFILE, fnm); + out_xvgrhistfile = opt2fn("-rhist", NFILE, fnm); + out_xvgkhistfile = opt2fn("-khist", NFILE, fnm); + out_xvginstefffile = opt2fn("-oe", NFILE, fnm); + out_datfile = opt2fn("-o",NFILE,fnm); + + bRKout = opt2bSet("-ot", NFILE, fnm); + bRhistout = opt2bSet("-rhist", NFILE, fnm); + bKhistout = opt2bSet("-khist", NFILE, fnm); + bDatout = opt2bSet("-o", NFILE, fnm); + bInstEffout = opt2bSet("-oe", NFILE, fnm); + + + /* PBC warning. */ + if (bPBCdist) + { + printf("Calculating distances to periodic image.\n"); + printf("Be careful! This produces only valid results for PBC in all three dimensions\n"); + } + + + if (bInstEffout && R0<=0.) + { + gmx_fatal(FARGS,"You have to specify R0 and R0 has to be larger than 0 nm.\n\n"); + } + + printf("Select group with donor atom pairs defining the transition moment\n"); + get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &ndon, &donindex,&grpnm); + + printf("Select group with acceptor atom pairs defining the transition moment\n"); + get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &nacc, &accindex,&grpnm); + + printf("Reading first frame\n"); + /* open trx file for reading */ + flags=0; + flags = flags | TRX_READ_X; + bHaveFirstFrame = read_first_frame(oenv, &status, in_trajfile, &fr, flags); + + if (bHaveFirstFrame) + { + printf("First frame is OK\n"); + natoms = fr.natoms; + if ((ndon % 2 != 0) || (nacc % 2 != 0)) + { + indexOK = FALSE; + } + else + { + for (i = 0; i < ndon;i++) + { + if (donindex[i] >= natoms) + indexOK = FALSE; + } + for (i = 0; i < nacc;i++) + { + if (accindex[i] >= natoms) + indexOK = FALSE; + } + } + + if (indexOK) + { + + if (bDatout) + { + datfp = fopen(out_datfile,"w"); + } + + if (bRKout) + { + rkfp = xvgropen(out_xvgrkfile, + "Distance and \\f{Symbol}k\\f{}\\S2\\N trajectory", + "Time (ps)", "Distance (nm) / \\f{Symbol}k\\f{}\\S2\\N", + oenv); + xvgr_legend(rkfp, 2, rkleg, oenv); + } + + if (bInstEffout) + { + iefp = xvgropen(out_xvginstefffile, + "Instantaneous RET Efficiency", + "Time (ps)", "RET Efficiency", + oenv); + xvgr_legend(iefp, 1, ieleg, oenv); + } + + + if (bRhistout) + { + snew(rvalues, allocblock); + rblocksallocated += 1; + snew(rhist, histbins); + } + + if (bKhistout) + { + snew(kappa2values, allocblock); + kblocksallocated += 1; + snew(khist, histbins); + } + + do + { + clear_rvec(donvec); + clear_rvec(accvec); + clear_rvec(donpos); + clear_rvec(accpos); + for (i = 0; i < ndon / 2; i++) + { + rvec_sub(donvec, fr.x[donindex[2 * i]], donvec); + rvec_add(donvec, fr.x[donindex[2 * i + 1]], donvec); + rvec_add(donpos, fr.x[donindex[2 * i]], donpos); + rvec_add(donpos, fr.x[donindex[2 * i + 1]], donpos); + } + + for (i = 0; i < nacc / 2; i++) + { + rvec_sub(accvec, fr.x[accindex[2 * i]], accvec); + rvec_add(accvec, fr.x[accindex[2 * i + 1]], accvec); + rvec_add(accpos, fr.x[accindex[2 * i]], accpos); + rvec_add(accpos, fr.x[accindex[2 * i + 1]], accpos); + } + + unitv(donvec, donvec); + unitv(accvec, accvec); + + svmul((real) 1. / ndon, donpos, donpos); + svmul((real) 1. / nacc, accpos, accpos); + + if (bPBCdist) + { + set_pbc(pbc, ePBC, fr.box); + pbc_dx(pbc, donpos, accpos, dist); + } + else + { + rvec_sub(donpos, accpos, dist); + } + + unitv(dist, distnorm); + R = norm(dist); + kappa2 = iprod(donvec, accvec)- 3.* (iprod(donvec, distnorm) * iprod(distnorm, accvec)); + kappa2 *= kappa2; + if (R0>0) + { + Rfrac=R/R0; + insteff=1/(1+(Rfrac*Rfrac*Rfrac*Rfrac*Rfrac*Rfrac)*2/3/kappa2); + insteffs+=insteff; + + if (bInstEffout) + { + fprintf(iefp, "%12.7f %12.7f\n", fr.time, insteff); + } + } + + + Rs += R; + kappa2s += kappa2; + rkcount++; + + if (bRKout) + fprintf(rkfp, "%12.7f %12.7f %12.7f\n", fr.time, R, kappa2); + + if (bDatout) + fprintf(datfp, "%12.7f %12.7f %12.7f\n", fr.time, R, kappa2); + + if (bRhistout) + { + rvalues[rkcount-1] = R; + if (rkcount % allocblock == 0) + { + srenew(rvalues, allocblock*(rblocksallocated+1)); + rblocksallocated += 1; + } + } + + if (bKhistout) + { + kappa2values[rkcount-1] = kappa2; + if (rkcount % allocblock == 0) + { + srenew(kappa2values, allocblock*(kblocksallocated+1)); + kblocksallocated += 1; + } + } + + bHaveNextFrame = read_next_frame(oenv, status, &fr); + } while (bHaveNextFrame); + + if (bRKout) + ffclose(rkfp); + + if (bDatout) + ffclose(datfp); + + if (bInstEffout) + ffclose(iefp); + + + if (bRhistout) + { + printf("Writing R-Histogram\n"); + rmin = rvalues[0]; + rmax = rvalues[0]; + for (i = 1; i < rkcount; i++) + { + if (rvalues[i] < rmin) + rmin = rvalues[i]; + else if (rvalues[i] > rmax) + rmax = rvalues[i]; + } + rmin -= histexpand; + rmax += histexpand; + + rrange = rmax - rmin; + rincr = rrange / histbins; + + for (i = 1; i < rkcount; i++) + { + bin = (int) ((rvalues[i] - rmin) / rincr); + rhist[bin] += 1; + } + if (bNormHist) + { + for (i = 0; i < histbins; i++) + rhist[i] /= rkcount * rrange/histbins; + rhfp = xvgropen(out_xvgrhistfile, "Distance Distribution", + "R (nm)", "Normalized Probability", oenv); + } else + { + rhfp = xvgropen(out_xvgrhistfile, "Distance Distribution", + "R (nm)", "Probability", oenv); + } + xvgr_legend(rhfp, 1, rhleg, oenv); + for (i = 0; i < histbins; i++) + { + fprintf(rhfp, "%12.7f %12.7f\n", (i + 0.5) * rincr + rmin, + rhist[i]); + } + ffclose(rhfp); + } + + if (bKhistout) + { + printf("Writing kappa^2-Histogram\n"); + krange = kmax - kmin; + kincr = krange / histbins; + + for (i = 1; i < rkcount; i++) + { + bin = (int) ((kappa2values[i] - kmin) / kincr); + khist[bin] += 1; + } + if (bNormHist) + { + for (i = 0; i < histbins; i++) + khist[i] /= rkcount * krange/histbins; + khfp = xvgropen(out_xvgkhistfile, + "\\f{Symbol}k\\f{}\\S2\\N Distribution", + "\\f{Symbol}k\\f{}\\S2\\N", + "Normalized Probability", oenv); + } else + { + khfp = xvgropen(out_xvgkhistfile, + "\\f{Symbol}k\\f{}\\S2\\N Distribution", + "\\f{Symbol}k\\f{}\\S2\\N", "Probability", oenv); + } + xvgr_legend(khfp, 1, khleg, oenv); + for (i = 0; i < histbins; i++) + { + fprintf(khfp, "%12.7f %12.7f\n", (i + 0.5) * kincr + kmin, + khist[i]); + } + ffclose(khfp); + } + + printf("\nAverages:\n"); + printf("R_avg = %8.4f nm\nKappa^2 = %8.4f\n", Rs / rkcount, + kappa2s / rkcount); + if (R0>0) + { + printf("E_RETavg = %8.4f\n", insteffs / rkcount); + } + please_cite(stdout,"Hoefling2011"); + } + else + { + gmx_fatal(FARGS,"Index file invalid, check your index file for correct pairs.\n"); + } + } + else + { + gmx_fatal(FARGS,"Could not read first frame of the trajectory.\n"); + } + + thanx(stderr); + return 0; +} + -- 2.22.0