src/tools/gmx_dyecoupl.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   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.
  13
  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.
  18  *
  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.
  25  *
  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.
  28  *
  29  * For more info, check our website at http://www.gromacs.org
  30  */
  31 #include <copyrite.h>
  32 #include <filenm.h>
  33 #include <macros.h>
  34 #include <pbc.h>
  35 #include <smalloc.h>
  36 #include <statutil.h>
  37 #include <vec.h>
  38 #include <xvgr.h>
  39 #include <nbsearch.h>
  40 #include <trajana.h>
  41 #include <math.h>
  42 #include "gmx_ana.h"
  43
  44
  45 int gmx_dyecoupl(int argc, char *argv[])
  46 {
  47     const char *desc[] =
  48     {
  49             "This tool extracts dye dynamics from trajectory files.",
  50             "Currently, R and kappa^2 between dyes is extracted for (F)RET",
  51             "simulations with assumed dipolar coupling as in the Foerster equation.",
  52             "It further allows the calculation of R(t) and kappa^2(t), R and",
  53             "kappa^2 histograms and averages, as well as the instantaneous FRET",
  54             "efficiency E(t) for a specified Foerster radius R_0 (switch [TT]-R0[tt]).",
  55             "The input dyes have to be whole (see res and mol pbc options",
  56             "in [TT]trjconv[tt]).",
  57             "The dye transition dipole moment has to be defined by at least",
  58             "a single atom pair, however multiple atom pairs can be provided ",
  59             "in the index file. The distance R is calculated on the basis of",
  60             "the COMs of the given atom pairs.",
  61             "The [TT]-pbcdist[tt] option calculates distances to the nearest periodic",
  62             "image instead to the distance in the box. This works however only,"
  63             "for periodic boundaries in all 3 dimensions.",
  64             "The [TT]-norm[tt] option (area-) normalizes the histograms."
  65     };
  66
  67         static gmx_bool bPBCdist = FALSE, bNormHist = FALSE;
  68     int histbins = 50;
  69     output_env_t oenv;
  70     real R0=-1;
  71
  72     t_pargs pa[] =
  73     {
  74             { "-pbcdist", FALSE, etBOOL, { &bPBCdist },"Distance R based on PBC" },
  75             { "-norm", FALSE, etBOOL, { &bNormHist },"Normalize histograms" },
  76             { "-bins", FALSE, etINT, {&histbins},"# of histogram bins" },
  77             { "-R0", FALSE, etREAL, {&R0},"Foerster radius including kappa^2=2/3 in nm" }
  78     };
  79 #define NPA asize(pa)
  80
  81     t_filenm fnm[] =
  82     {
  83             { efTRX, "-f", NULL, ffREAD },
  84             { efNDX, NULL, NULL, ffREAD },
  85             { efXVG, "-ot", "rkappa",ffOPTWR },
  86             { efXVG, "-oe", "insteff",ffOPTWR },
  87             { efDAT, "-o", "rkappa",ffOPTWR },
  88             { efXVG, "-rhist","rhist", ffOPTWR },
  89             { efXVG, "-khist", "khist", ffOPTWR }
  90     };
  91 #define NFILE asize(fnm)
  92
  93
  94     const char *in_trajfile, *in_ndxfile, *out_xvgrkfile = NULL, *out_xvginstefffile = NULL, *out_xvgrhistfile = NULL, *out_xvgkhistfile = NULL,*out_datfile=NULL;
  95     gmx_bool bHaveFirstFrame, bHaveNextFrame, indexOK = TRUE;
  96     int ndon, nacc;
  97     atom_id *donindex, *accindex;
  98     char *grpnm;
  99     t_atoms *atoms = NULL;
 100     t_trxstatus *status;
 101     t_trxframe fr;
 102
 103     int flags;
 104     int allocblock = 1000;
 105     real histexpand = 1e-6;
 106     rvec donvec, accvec, donpos, accpos, dist, distnorm;
 107     int natoms;
 108
 109     /*we rely on PBC autodetection (...currently)*/
 110     int ePBC = -1;
 111
 112     real *rvalues=NULL, *kappa2values=NULL, *rhist=NULL, *khist=NULL;
 113     t_pbc *pbc=NULL;
 114     int i, bin;
 115     FILE *rkfp = NULL, *rhfp = NULL, *khfp = NULL,*datfp=NULL,*iefp=NULL;
 116     gmx_bool bRKout, bRhistout, bKhistout,bDatout,bInstEffout,grident;
 117
 118     const char *rkleg[2] = { "R", "\\f{Symbol}k\\f{}\\S2\\N" };
 119     const char *rhleg[1] = { "p(R)" };
 120     const char *khleg[1] = { "p(\\f{Symbol}k\\f{}\\S2\\N)" };
 121     const char *ieleg[1] = { "E\\sRET\\N(t)" };
 122
 123     real R, kappa2, insteff, Rs = 0., kappa2s = 0., insteffs=0., rmax, rmin, kmin = 0., kmax = 4.,
 124             rrange, krange, rincr, kincr,Rfrac;
 125     int rkcount = 0, rblocksallocated = 0, kblocksallocated = 0;
 126
 127     CopyRight(stderr, argv[0]);
 128     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
 130
 131     /* Check command line options for filenames and set bool flags when switch used*/
 132     in_trajfile = opt2fn("-f", NFILE, fnm);
 133     in_ndxfile = opt2fn("-n", NFILE, fnm);
 134     out_xvgrkfile = opt2fn("-ot", NFILE, fnm);
 135     out_xvgrhistfile = opt2fn("-rhist", NFILE, fnm);
 136     out_xvgkhistfile = opt2fn("-khist", NFILE, fnm);
 137     out_xvginstefffile = opt2fn("-oe", NFILE, fnm);
 138     out_datfile = opt2fn("-o",NFILE,fnm);
 139
 140     bRKout = opt2bSet("-ot", NFILE, fnm);
 141     bRhistout = opt2bSet("-rhist", NFILE, fnm);
 142     bKhistout = opt2bSet("-khist", NFILE, fnm);
 143     bDatout = opt2bSet("-o", NFILE, fnm);
 144     bInstEffout = opt2bSet("-oe", NFILE, fnm);
 145
 146
 147     /* PBC warning. */
 148     if (bPBCdist)
 149     {
 150         printf("Calculating distances to periodic image.\n");
 151         printf("Be careful! This produces only valid results for PBC in all three dimensions\n");
 152     }
 153
 154
 155     if (bInstEffout && R0<=0.)
 156     {
 157         gmx_fatal(FARGS,"You have to specify R0 and R0 has to be larger than 0 nm.\n\n");
 158     }
 159
 160     printf("Select group with donor atom pairs defining the transition moment\n");
 161     get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &ndon, &donindex,&grpnm);
 162
 163     printf("Select group with acceptor atom pairs defining the transition moment\n");
 164     get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &nacc, &accindex,&grpnm);
 165
 166     /*check if groups are identical*/
 167     grident=TRUE;
 168
 169     if (ndon == nacc)
 170     {
 171         for (i=0;i<nacc;i++)
 172         {
 173             if(accindex[i] != donindex[i])
 174             {
 175                 grident=FALSE;
 176                 break;
 177             }
 178         }
 179     }
 180
 181     if (grident)
 182     {
 183         gmx_fatal(FARGS,"Donor and acceptor group are identical. This makes no sense.");
 184     }
 185
 186     printf("Reading first frame\n");
 187     /* open trx file for reading */
 188     flags=0;
 189     flags = flags | TRX_READ_X;
 190     bHaveFirstFrame = read_first_frame(oenv, &status, in_trajfile, &fr, flags);
 191
 192     if (bHaveFirstFrame)
 193     {
 194         printf("First frame is OK\n");
 195         natoms = fr.natoms;
 196         if ((ndon % 2 != 0) || (nacc % 2 != 0))
 197         {
 198             indexOK = FALSE;
 199         }
 200         else
 201         {
 202             for (i = 0; i < ndon;i++)
 203             {
 204                 if (donindex[i] >= natoms)
 205                     indexOK = FALSE;
 206             }
 207             for (i = 0; i < nacc;i++)
 208             {
 209                 if (accindex[i] >= natoms)
 210                     indexOK = FALSE;
 211             }
 212         }
 213
 214         if (indexOK)
 215         {
 216
 217             if (bDatout)
 218             {
 219                 datfp = fopen(out_datfile,"w");
 220             }
 221
 222             if (bRKout)
 223             {
 224                 rkfp = xvgropen(out_xvgrkfile,
 225                         "Distance and \\f{Symbol}k\\f{}\\S2\\N trajectory",
 226                         "Time (ps)", "Distance (nm) / \\f{Symbol}k\\f{}\\S2\\N",
 227                         oenv);
 228                 xvgr_legend(rkfp, 2, rkleg, oenv);
 229             }
 230
 231             if (bInstEffout)
 232             {
 233                 iefp = xvgropen(out_xvginstefffile,
 234                         "Instantaneous RET Efficiency",
 235                         "Time (ps)", "RET Efficiency",
 236                         oenv);
 237                 xvgr_legend(iefp, 1, ieleg, oenv);
 238             }
 239
 240
 241             if (bRhistout)
 242             {
 243                 snew(rvalues, allocblock);
 244                 rblocksallocated += 1;
 245                 snew(rhist, histbins);
 246             }
 247
 248             if (bKhistout)
 249             {
 250                 snew(kappa2values, allocblock);
 251                 kblocksallocated += 1;
 252                 snew(khist, histbins);
 253             }
 254
 255             do
 256             {
 257                 clear_rvec(donvec);
 258                 clear_rvec(accvec);
 259                 clear_rvec(donpos);
 260                 clear_rvec(accpos);
 261                 for (i = 0; i < ndon / 2; i++)
 262                 {
 263                     rvec_sub(donvec, fr.x[donindex[2 * i]], donvec);
 264                     rvec_add(donvec, fr.x[donindex[2 * i + 1]], donvec);
 265                     rvec_add(donpos, fr.x[donindex[2 * i]], donpos);
 266                     rvec_add(donpos, fr.x[donindex[2 * i + 1]], donpos);
 267                 }
 268
 269                 for (i = 0; i < nacc / 2; i++)
 270                 {
 271                     rvec_sub(accvec, fr.x[accindex[2 * i]], accvec);
 272                     rvec_add(accvec, fr.x[accindex[2 * i + 1]], accvec);
 273                     rvec_add(accpos, fr.x[accindex[2 * i]], accpos);
 274                     rvec_add(accpos, fr.x[accindex[2 * i + 1]], accpos);
 275                 }
 276
 277                 unitv(donvec, donvec);
 278                 unitv(accvec, accvec);
 279
 280                 svmul((real) 1. / ndon, donpos, donpos);
 281                 svmul((real) 1. / nacc, accpos, accpos);
 282
 283                 if (bPBCdist)
 284                 {
 285                     set_pbc(pbc, ePBC, fr.box);
 286                     pbc_dx(pbc, donpos, accpos, dist);
 287                 }
 288                 else
 289                 {
 290                     rvec_sub(donpos, accpos, dist);
 291                 }
 292
 293                 unitv(dist, distnorm);
 294                 R = norm(dist);
 295                 kappa2 = iprod(donvec, accvec)- 3.* (iprod(donvec, distnorm) * iprod(distnorm, accvec));
 296                 kappa2 *= kappa2;
 297                 if (R0>0)
 298                 {
 299                     Rfrac=R/R0;
 300                     insteff=1/(1+(Rfrac*Rfrac*Rfrac*Rfrac*Rfrac*Rfrac)*2/3/kappa2);
 301                     insteffs+=insteff;
 302
 303                     if (bInstEffout)
 304                     {
 305                         fprintf(iefp, "%12.7f %12.7f\n", fr.time, insteff);
 306                     }
 307                 }
 308
 309
 310                 Rs += R;
 311                 kappa2s += kappa2;
 312                 rkcount++;
 313
 314                 if (bRKout)
 315                     fprintf(rkfp, "%12.7f %12.7f %12.7f\n", fr.time, R, kappa2);
 316
 317                 if (bDatout)
 318                     fprintf(datfp, "%12.7f %12.7f %12.7f\n", fr.time, R, kappa2);
 319
 320                 if (bRhistout)
 321                 {
 322                     rvalues[rkcount-1] = R;
 323                     if (rkcount % allocblock == 0)
 324                     {
 325                         srenew(rvalues, allocblock*(rblocksallocated+1));
 326                         rblocksallocated += 1;
 327                     }
 328                 }
 329
 330                 if (bKhistout)
 331                 {
 332                     kappa2values[rkcount-1] = kappa2;
 333                     if (rkcount % allocblock == 0)
 334                     {
 335                         srenew(kappa2values, allocblock*(kblocksallocated+1));
 336                         kblocksallocated += 1;
 337                     }
 338                 }
 339
 340                 bHaveNextFrame = read_next_frame(oenv, status, &fr);
 341             } while (bHaveNextFrame);
 342
 343             if (bRKout)
 344                 ffclose(rkfp);
 345
 346             if (bDatout)
 347                 ffclose(datfp);
 348
 349             if (bInstEffout)
 350                 ffclose(iefp);
 351
 352
 353             if (bRhistout)
 354             {
 355                 printf("Writing R-Histogram\n");
 356                 rmin = rvalues[0];
 357                 rmax = rvalues[0];
 358                 for (i = 1; i < rkcount; i++)
 359                 {
 360                     if (rvalues[i] < rmin)
 361                         rmin = rvalues[i];
 362                     else if (rvalues[i] > rmax)
 363                         rmax = rvalues[i];
 364                 }
 365                 rmin -= histexpand;
 366                 rmax += histexpand;
 367
 368                 rrange = rmax - rmin;
 369                 rincr = rrange / histbins;
 370
 371                 for (i = 1; i < rkcount; i++)
 372                 {
 373                     bin = (int) ((rvalues[i] - rmin) / rincr);
 374                     rhist[bin] += 1;
 375                 }
 376                 if (bNormHist)
 377                 {
 378                     for (i = 0; i < histbins; i++)
 379                         rhist[i] /= rkcount * rrange/histbins;
 380                     rhfp = xvgropen(out_xvgrhistfile, "Distance Distribution",
 381                             "R (nm)", "Normalized Probability", oenv);
 382                 } else
 383                 {
 384                     rhfp = xvgropen(out_xvgrhistfile, "Distance Distribution",
 385                             "R (nm)", "Probability", oenv);
 386                 }
 387                 xvgr_legend(rhfp, 1, rhleg, oenv);
 388                 for (i = 0; i < histbins; i++)
 389                 {
 390                     fprintf(rhfp, "%12.7f %12.7f\n", (i + 0.5) * rincr + rmin,
 391                             rhist[i]);
 392                 }
 393                 ffclose(rhfp);
 394             }
 395
 396             if (bKhistout)
 397             {
 398                 printf("Writing kappa^2-Histogram\n");
 399                 krange = kmax - kmin;
 400                 kincr = krange / histbins;
 401
 402                 for (i = 1; i < rkcount; i++)
 403                 {
 404                     bin = (int) ((kappa2values[i] - kmin) / kincr);
 405                     khist[bin] += 1;
 406                 }
 407                 if (bNormHist)
 408                 {
 409                     for (i = 0; i < histbins; i++)
 410                         khist[i] /= rkcount * krange/histbins;
 411                     khfp = xvgropen(out_xvgkhistfile,
 412                             "\\f{Symbol}k\\f{}\\S2\\N Distribution",
 413                             "\\f{Symbol}k\\f{}\\S2\\N",
 414                             "Normalized Probability", oenv);
 415                 } else
 416                 {
 417                     khfp = xvgropen(out_xvgkhistfile,
 418                             "\\f{Symbol}k\\f{}\\S2\\N Distribution",
 419                             "\\f{Symbol}k\\f{}\\S2\\N", "Probability", oenv);
 420                 }
 421                 xvgr_legend(khfp, 1, khleg, oenv);
 422                 for (i = 0; i < histbins; i++)
 423                 {
 424                     fprintf(khfp, "%12.7f %12.7f\n", (i + 0.5) * kincr + kmin,
 425                             khist[i]);
 426                 }
 427                 ffclose(khfp);
 428             }
 429
 430             printf("\nAverages:\n");
 431             printf("R_avg   = %8.4f nm\nKappa^2 = %8.4f\n", Rs / rkcount,
 432                     kappa2s / rkcount);
 433             if (R0>0)
 434             {
 435                 printf("E_RETavg   = %8.4f\n", insteffs / rkcount);
 436             }
 437             please_cite(stdout,"Hoefling2011");
 438         }
 439         else
 440         {
 441             gmx_fatal(FARGS,"Index file invalid, check your index file for correct pairs.\n");
 442         }
 443     }
 444     else
 445     {
 446         gmx_fatal(FARGS,"Could not read first frame of the trajectory.\n");
 447     }
 448
 449     thanx(stderr);
 450     return 0;
 451 }
 452