src/tools/gmx_dyecoupl.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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, 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.
  11  *
  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.
  16  *
  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.
  21  *
  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.
  26  *
  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.
  34  *
  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.
  37  */
  38 #include <copyrite.h>
  39 #include <filenm.h>
  40 #include <macros.h>
  41 #include <pbc.h>
  42 #include <smalloc.h>
  43 #include <statutil.h>
  44 #include <vec.h>
  45 #include <xvgr.h>
  46 #include <nbsearch.h>
  47 #include <trajana.h>
  48 #include <math.h>
  49 #include "gmx_ana.h"
  50
  51
  52 int gmx_dyecoupl(int argc, char *argv[])
  53 {
  54     const char *desc[] =
  55     {
  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."
  72     };
  73
  74         static gmx_bool bPBCdist = FALSE, bNormHist = FALSE;
  75     int histbins = 50;
  76     output_env_t oenv;
  77     real R0=-1;
  78
  79     t_pargs pa[] =
  80     {
  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" }
  85     };
  86 #define NPA asize(pa)
  87
  88     t_filenm fnm[] =
  89     {
  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 }
  97     };
  98 #define NFILE asize(fnm)
  99
 100
 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;
 103     int ndon, nacc;
 104     atom_id *donindex, *accindex;
 105     char *grpnm;
 106     t_atoms *atoms = NULL;
 107     t_trxstatus *status;
 108     t_trxframe fr;
 109
 110     int flags;
 111     int allocblock = 1000;
 112     real histexpand = 1e-6;
 113     rvec donvec, accvec, donpos, accpos, dist, distnorm;
 114     int natoms;
 115
 116     /*we rely on PBC autodetection (...currently)*/
 117     int ePBC = -1;
 118
 119     real *rvalues=NULL, *kappa2values=NULL, *rhist=NULL, *khist=NULL;
 120     t_pbc *pbc=NULL;
 121     int i, bin;
 122     FILE *rkfp = NULL, *rhfp = NULL, *khfp = NULL,*datfp=NULL,*iefp=NULL;
 123     gmx_bool bRKout, bRhistout, bKhistout,bDatout,bInstEffout,grident;
 124
 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)" };
 129
 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;
 133
 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);
 136
 137
 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);
 146
 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);
 152
 153
 154     /* PBC warning. */
 155     if (bPBCdist)
 156     {
 157         printf("Calculating distances to periodic image.\n");
 158         printf("Be careful! This produces only valid results for PBC in all three dimensions\n");
 159     }
 160
 161
 162     if (bInstEffout && R0<=0.)
 163     {
 164         gmx_fatal(FARGS,"You have to specify R0 and R0 has to be larger than 0 nm.\n\n");
 165     }
 166
 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);
 169
 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);
 172
 173     /*check if groups are identical*/
 174     grident=TRUE;
 175
 176     if (ndon == nacc)
 177     {
 178         for (i=0;i<nacc;i++)
 179         {
 180             if(accindex[i] != donindex[i])
 181             {
 182                 grident=FALSE;
 183                 break;
 184             }
 185         }
 186     }
 187
 188     if (grident)
 189     {
 190         gmx_fatal(FARGS,"Donor and acceptor group are identical. This makes no sense.");
 191     }
 192
 193     printf("Reading first frame\n");
 194     /* open trx file for reading */
 195     flags=0;
 196     flags = flags | TRX_READ_X;
 197     bHaveFirstFrame = read_first_frame(oenv, &status, in_trajfile, &fr, flags);
 198
 199     if (bHaveFirstFrame)
 200     {
 201         printf("First frame is OK\n");
 202         natoms = fr.natoms;
 203         if ((ndon % 2 != 0) || (nacc % 2 != 0))
 204         {
 205             indexOK = FALSE;
 206         }
 207         else
 208         {
 209             for (i = 0; i < ndon;i++)
 210             {
 211                 if (donindex[i] >= natoms)
 212                     indexOK = FALSE;
 213             }
 214             for (i = 0; i < nacc;i++)
 215             {
 216                 if (accindex[i] >= natoms)
 217                     indexOK = FALSE;
 218             }
 219         }
 220
 221         if (indexOK)
 222         {
 223
 224             if (bDatout)
 225             {
 226                 datfp = fopen(out_datfile,"w");
 227             }
 228
 229             if (bRKout)
 230             {
 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",
 234                         oenv);
 235                 xvgr_legend(rkfp, 2, rkleg, oenv);
 236             }
 237
 238             if (bInstEffout)
 239             {
 240                 iefp = xvgropen(out_xvginstefffile,
 241                         "Instantaneous RET Efficiency",
 242                         "Time (ps)", "RET Efficiency",
 243                         oenv);
 244                 xvgr_legend(iefp, 1, ieleg, oenv);
 245             }
 246
 247
 248             if (bRhistout)
 249             {
 250                 snew(rvalues, allocblock);
 251                 rblocksallocated += 1;
 252                 snew(rhist, histbins);
 253             }
 254
 255             if (bKhistout)
 256             {
 257                 snew(kappa2values, allocblock);
 258                 kblocksallocated += 1;
 259                 snew(khist, histbins);
 260             }
 261
 262             do
 263             {
 264                 clear_rvec(donvec);
 265                 clear_rvec(accvec);
 266                 clear_rvec(donpos);
 267                 clear_rvec(accpos);
 268                 for (i = 0; i < ndon / 2; i++)
 269                 {
 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);
 274                 }
 275
 276                 for (i = 0; i < nacc / 2; i++)
 277                 {
 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);
 282                 }
 283
 284                 unitv(donvec, donvec);
 285                 unitv(accvec, accvec);
 286
 287                 svmul((real) 1. / ndon, donpos, donpos);
 288                 svmul((real) 1. / nacc, accpos, accpos);
 289
 290                 if (bPBCdist)
 291                 {
 292                     set_pbc(pbc, ePBC, fr.box);
 293                     pbc_dx(pbc, donpos, accpos, dist);
 294                 }
 295                 else
 296                 {
 297                     rvec_sub(donpos, accpos, dist);
 298                 }
 299
 300                 unitv(dist, distnorm);
 301                 R = norm(dist);
 302                 kappa2 = iprod(donvec, accvec)- 3.* (iprod(donvec, distnorm) * iprod(distnorm, accvec));
 303                 kappa2 *= kappa2;
 304                 if (R0>0)
 305                 {
 306                     Rfrac=R/R0;
 307                     insteff=1/(1+(Rfrac*Rfrac*Rfrac*Rfrac*Rfrac*Rfrac)*2/3/kappa2);
 308                     insteffs+=insteff;
 309
 310                     if (bInstEffout)
 311                     {
 312                         fprintf(iefp, "%12.7f %12.7f\n", fr.time, insteff);
 313                     }
 314                 }
 315
 316
 317                 Rs += R;
 318                 kappa2s += kappa2;
 319                 rkcount++;
 320
 321                 if (bRKout)
 322                     fprintf(rkfp, "%12.7f %12.7f %12.7f\n", fr.time, R, kappa2);
 323
 324                 if (bDatout)
 325                     fprintf(datfp, "%12.7f %12.7f %12.7f\n", fr.time, R, kappa2);
 326
 327                 if (bRhistout)
 328                 {
 329                     rvalues[rkcount-1] = R;
 330                     if (rkcount % allocblock == 0)
 331                     {
 332                         srenew(rvalues, allocblock*(rblocksallocated+1));
 333                         rblocksallocated += 1;
 334                     }
 335                 }
 336
 337                 if (bKhistout)
 338                 {
 339                     kappa2values[rkcount-1] = kappa2;
 340                     if (rkcount % allocblock == 0)
 341                     {
 342                         srenew(kappa2values, allocblock*(kblocksallocated+1));
 343                         kblocksallocated += 1;
 344                     }
 345                 }
 346
 347                 bHaveNextFrame = read_next_frame(oenv, status, &fr);
 348             } while (bHaveNextFrame);
 349
 350             if (bRKout)
 351                 ffclose(rkfp);
 352
 353             if (bDatout)
 354                 ffclose(datfp);
 355
 356             if (bInstEffout)
 357                 ffclose(iefp);
 358
 359
 360             if (bRhistout)
 361             {
 362                 printf("Writing R-Histogram\n");
 363                 rmin = rvalues[0];
 364                 rmax = rvalues[0];
 365                 for (i = 1; i < rkcount; i++)
 366                 {
 367                     if (rvalues[i] < rmin)
 368                         rmin = rvalues[i];
 369                     else if (rvalues[i] > rmax)
 370                         rmax = rvalues[i];
 371                 }
 372                 rmin -= histexpand;
 373                 rmax += histexpand;
 374
 375                 rrange = rmax - rmin;
 376                 rincr = rrange / histbins;
 377
 378                 for (i = 1; i < rkcount; i++)
 379                 {
 380                     bin = (int) ((rvalues[i] - rmin) / rincr);
 381                     rhist[bin] += 1;
 382                 }
 383                 if (bNormHist)
 384                 {
 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);
 389                 } else
 390                 {
 391                     rhfp = xvgropen(out_xvgrhistfile, "Distance Distribution",
 392                             "R (nm)", "Probability", oenv);
 393                 }
 394                 xvgr_legend(rhfp, 1, rhleg, oenv);
 395                 for (i = 0; i < histbins; i++)
 396                 {
 397                     fprintf(rhfp, "%12.7f %12.7f\n", (i + 0.5) * rincr + rmin,
 398                             rhist[i]);
 399                 }
 400                 ffclose(rhfp);
 401             }
 402
 403             if (bKhistout)
 404             {
 405                 printf("Writing kappa^2-Histogram\n");
 406                 krange = kmax - kmin;
 407                 kincr = krange / histbins;
 408
 409                 for (i = 1; i < rkcount; i++)
 410                 {
 411                     bin = (int) ((kappa2values[i] - kmin) / kincr);
 412                     khist[bin] += 1;
 413                 }
 414                 if (bNormHist)
 415                 {
 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);
 422                 } else
 423                 {
 424                     khfp = xvgropen(out_xvgkhistfile,
 425                             "\\f{Symbol}k\\f{}\\S2\\N Distribution",
 426                             "\\f{Symbol}k\\f{}\\S2\\N", "Probability", oenv);
 427                 }
 428                 xvgr_legend(khfp, 1, khleg, oenv);
 429                 for (i = 0; i < histbins; i++)
 430                 {
 431                     fprintf(khfp, "%12.7f %12.7f\n", (i + 0.5) * kincr + kmin,
 432                             khist[i]);
 433                 }
 434                 ffclose(khfp);
 435             }
 436
 437             printf("\nAverages:\n");
 438             printf("R_avg   = %8.4f nm\nKappa^2 = %8.4f\n", Rs / rkcount,
 439                     kappa2s / rkcount);
 440             if (R0>0)
 441             {
 442                 printf("E_RETavg   = %8.4f\n", insteffs / rkcount);
 443             }
 444             please_cite(stdout,"Hoefling2011");
 445         }
 446         else
 447         {
 448             gmx_fatal(FARGS,"Index file invalid, check your index file for correct pairs.\n");
 449         }
 450     }
 451     else
 452     {
 453         gmx_fatal(FARGS,"Could not read first frame of the trajectory.\n");
 454     }
 455
 456     thanx(stderr);
 457     return 0;
 458 }
 459