src/gromacs/gmxana/gmx_kinetics.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  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38 #include <math.h>
  39 #include <string.h>
  40 #include <stdlib.h>
  41 #include "statutil.h"
  42 #include "sysstuff.h"
  43 #include "typedefs.h"
  44 #include "smalloc.h"
  45 #include "macros.h"
  46 #include "gmx_fatal.h"
  47 #include "vec.h"
  48 #include "copyrite.h"
  49 #include "futil.h"
  50 #include "readinp.h"
  51 #include "statutil.h"
  52 #include "txtdump.h"
  53 #include "gstat.h"
  54 #include "xvgr.h"
  55 #include "physics.h"
  56 #include "gmx_ana.h"
  57
  58 enum {
  59     epAuf, epEuf, epAfu, epEfu, epNR
  60 };
  61 enum {
  62     eqAif, eqEif, eqAfi, eqEfi, eqAui, eqEui, eqAiu, eqEiu, eqNR
  63 };
  64 static char *eep[epNR] = { "Af", "Ef", "Au", "Eu" };
  65 static char *eeq[eqNR] = { "Aif", "Eif", "Afi", "Efi", "Aui", "Eui", "Aiu", "Eiu" };
  66
  67 typedef struct {
  68     int       nreplica;  /* Number of replicas in the calculation                   */
  69     int       nframe;    /* Number of time frames                                   */
  70     int       nstate;    /* Number of states the system can be in, e.g. F,I,U       */
  71     int       nparams;   /* Is 2, 4 or 8                                            */
  72     gmx_bool *bMask;     /* Determine whether this replica is part of the d2 comp.  */
  73     gmx_bool  bSum;
  74     gmx_bool  bDiscrete; /* Use either discrete folding (0/1) or a continuous       */
  75     /* criterion */
  76     int       nmask;     /* Number of replicas taken into account                   */
  77     real      dt;        /* Timestep between frames                                 */
  78     int       j0, j1;    /* Range of frames used in calculating delta               */
  79     real    **temp, **data, **data2;
  80     int     **state;     /* State index running from 0 (F) to nstate-1 (U)          */
  81     real    **beta, **fcalt, **icalt;
  82     real     *time, *sumft, *sumit, *sumfct, *sumict;
  83     real     *params;
  84     real     *d2_replica;
  85 } t_remd_data;
  86
  87 #ifdef HAVE_LIBGSL
  88 #include <gsl/gsl_multimin.h>
  89
  90 static char *itoa(int i)
  91 {
  92     static char ptr[12];
  93
  94     sprintf(ptr, "%d", i);
  95     return ptr;
  96 }
  97
  98 static char *epnm(int nparams, int index)
  99 {
 100     static char buf[32], from[8], to[8];
 101     int         nn, ni, ii;
 102
 103     range_check(index, 0, nparams);
 104     if ((nparams == 2) || (nparams == 4))
 105     {
 106         return eep[index];
 107     }
 108     else if ((nparams > 4) && (nparams % 4 == 0))
 109     {
 110         return eeq[index];
 111     }
 112     else
 113     {
 114         gmx_fatal(FARGS, "Don't know how to handle %d parameters", nparams);
 115     }
 116
 117     return NULL;
 118 }
 119
 120 static gmx_bool bBack(t_remd_data *d)
 121 {
 122     return (d->nparams > 2);
 123 }
 124
 125 static real is_folded(t_remd_data *d, int irep, int jframe)
 126 {
 127     if (d->state[irep][jframe] == 0)
 128     {
 129         return 1.0;
 130     }
 131     else
 132     {
 133         return 0.0;
 134     }
 135 }
 136
 137 static real is_unfolded(t_remd_data *d, int irep, int jframe)
 138 {
 139     if (d->state[irep][jframe] == d->nstate-1)
 140     {
 141         return 1.0;
 142     }
 143     else
 144     {
 145         return 0.0;
 146     }
 147 }
 148
 149 static real is_intermediate(t_remd_data *d, int irep, int jframe)
 150 {
 151     if ((d->state[irep][jframe] == 1) && (d->nstate > 2))
 152     {
 153         return 1.0;
 154     }
 155     else
 156     {
 157         return 0.0;
 158     }
 159 }
 160
 161 static void integrate_dfdt(t_remd_data *d)
 162 {
 163     int    i, j;
 164     double beta, ddf, ddi, df, db, fac, sumf, sumi, area;
 165
 166     d->sumfct[0] = 0;
 167     d->sumict[0] = 0;
 168     for (i = 0; (i < d->nreplica); i++)
 169     {
 170         if (d->bMask[i])
 171         {
 172             if (d->bDiscrete)
 173             {
 174                 ddf = 0.5*d->dt*is_folded(d, i, 0);
 175                 ddi = 0.5*d->dt*is_intermediate(d, i, 0);
 176             }
 177             else
 178             {
 179                 ddf = 0.5*d->dt*d->state[i][0];
 180                 ddi = 0.0;
 181             }
 182             d->fcalt[i][0] = ddf;
 183             d->icalt[i][0] = ddi;
 184             d->sumfct[0]  += ddf;
 185             d->sumict[0]  += ddi;
 186         }
 187     }
 188     for (j = 1; (j < d->nframe); j++)
 189     {
 190         if (j == d->nframe-1)
 191         {
 192             fac = 0.5*d->dt;
 193         }
 194         else
 195         {
 196             fac = d->dt;
 197         }
 198         sumf = sumi = 0;
 199         for (i = 0; (i < d->nreplica); i++)
 200         {
 201             if (d->bMask[i])
 202             {
 203                 beta = d->beta[i][j];
 204                 if ((d->nstate <= 2) || d->bDiscrete)
 205                 {
 206                     if (d->bDiscrete)
 207                     {
 208                         df = (d->params[epAuf]*exp(-beta*d->params[epEuf])*
 209                               is_unfolded(d, i, j));
 210                     }
 211                     else
 212                     {
 213                         area = (d->data2 ? d->data2[i][j] : 1.0);
 214                         df   =  area*d->params[epAuf]*exp(-beta*d->params[epEuf]);
 215                     }
 216                     if (bBack(d))
 217                     {
 218                         db = 0;
 219                         if (d->bDiscrete)
 220                         {
 221                             db = (d->params[epAfu]*exp(-beta*d->params[epEfu])*
 222                                   is_folded(d, i, j));
 223                         }
 224                         else
 225                         {
 226                             gmx_fatal(FARGS, "Back reaction not implemented with continuous");
 227                         }
 228                         ddf = fac*(df-db);
 229                     }
 230                     else
 231                     {
 232                         ddf = fac*df;
 233                     }
 234                     d->fcalt[i][j] = d->fcalt[i][j-1] + ddf;
 235                     sumf          += ddf;
 236                 }
 237                 else
 238                 {
 239                     ddf = fac*((d->params[eqAif]*exp(-beta*d->params[eqEif])*
 240                                 is_intermediate(d, i, j)) -
 241                                (d->params[eqAfi]*exp(-beta*d->params[eqEfi])*
 242                                 is_folded(d, i, j)));
 243                     ddi = fac*((d->params[eqAui]*exp(-beta*d->params[eqEui])*
 244                                 is_unfolded(d, i, j)) -
 245                                (d->params[eqAiu]*exp(-beta*d->params[eqEiu])*
 246                                 is_intermediate(d, i, j)));
 247                     d->fcalt[i][j] = d->fcalt[i][j-1] + ddf;
 248                     d->icalt[i][j] = d->icalt[i][j-1] + ddi;
 249                     sumf          += ddf;
 250                     sumi          += ddi;
 251                 }
 252             }
 253         }
 254         d->sumfct[j] = d->sumfct[j-1] + sumf;
 255         d->sumict[j] = d->sumict[j-1] + sumi;
 256     }
 257     if (debug)
 258     {
 259         fprintf(debug, "@type xy\n");
 260         for (j = 0; (j < d->nframe); j++)
 261         {
 262             fprintf(debug, "%8.3f  %12.5e\n", d->time[j], d->sumfct[j]);
 263         }
 264         fprintf(debug, "&\n");
 265     }
 266 }
 267
 268 static void sum_ft(t_remd_data *d)
 269 {
 270     int    i, j;
 271     double fac;
 272
 273     for (j = 0; (j < d->nframe); j++)
 274     {
 275         d->sumft[j] = 0;
 276         d->sumit[j] = 0;
 277         if ((j == 0) || (j == d->nframe-1))
 278         {
 279             fac = d->dt*0.5;
 280         }
 281         else
 282         {
 283             fac = d->dt;
 284         }
 285         for (i = 0; (i < d->nreplica); i++)
 286         {
 287             if (d->bMask[i])
 288             {
 289                 if (d->bDiscrete)
 290                 {
 291                     d->sumft[j] += fac*is_folded(d, i, j);
 292                     d->sumit[j] += fac*is_intermediate(d, i, j);
 293                 }
 294                 else
 295                 {
 296                     d->sumft[j] += fac*d->state[i][j];
 297                 }
 298             }
 299         }
 300     }
 301 }
 302
 303 static double calc_d2(t_remd_data *d)
 304 {
 305     int    i, j;
 306     double dd2, d2 = 0, dr2, tmp;
 307
 308     integrate_dfdt(d);
 309
 310     if (d->bSum)
 311     {
 312         for (j = d->j0; (j < d->j1); j++)
 313         {
 314             if (d->bDiscrete)
 315             {
 316                 d2  += sqr(d->sumft[j]-d->sumfct[j]);
 317                 if (d->nstate > 2)
 318                 {
 319                     d2 += sqr(d->sumit[j]-d->sumict[j]);
 320                 }
 321             }
 322             else
 323             {
 324                 d2  += sqr(d->sumft[j]-d->sumfct[j]);
 325             }
 326         }
 327     }
 328     else
 329     {
 330         for (i = 0; (i < d->nreplica); i++)
 331         {
 332             dr2 = 0;
 333             if (d->bMask[i])
 334             {
 335                 for (j = d->j0; (j < d->j1); j++)
 336                 {
 337                     tmp  = sqr(is_folded(d, i, j)-d->fcalt[i][j]);
 338                     d2  += tmp;
 339                     dr2 += tmp;
 340                     if (d->nstate > 2)
 341                     {
 342                         tmp  = sqr(is_intermediate(d, i, j)-d->icalt[i][j]);
 343                         d2  += tmp;
 344                         dr2 += tmp;
 345                     }
 346                 }
 347                 d->d2_replica[i] = dr2/(d->j1-d->j0);
 348             }
 349         }
 350     }
 351     dd2 = (d2/(d->j1-d->j0))/(d->bDiscrete ? d->nmask : 1);
 352
 353     return dd2;
 354 }
 355
 356 static double my_f(const gsl_vector *v, void *params)
 357 {
 358     t_remd_data *d       = (t_remd_data *) params;
 359     double       penalty = 0;
 360     int          i;
 361
 362     for (i = 0; (i < d->nparams); i++)
 363     {
 364         d->params[i] = gsl_vector_get(v, i);
 365         if (d->params[i] < 0)
 366         {
 367             penalty += 10;
 368         }
 369     }
 370     if (penalty > 0)
 371     {
 372         return penalty;
 373     }
 374     else
 375     {
 376         return calc_d2(d);
 377     }
 378 }
 379
 380 static void optimize_remd_parameters(t_remd_data *d, int maxiter,
 381                                      real tol)
 382 {
 383     real   size, d2;
 384     int    iter   = 0;
 385     int    status = 0;
 386     int    i;
 387
 388     const gsl_multimin_fminimizer_type *T;
 389     gsl_multimin_fminimizer            *s;
 390
 391     gsl_vector                         *x, *dx;
 392     gsl_multimin_function               my_func;
 393
 394     my_func.f      = &my_f;
 395     my_func.n      = d->nparams;
 396     my_func.params = (void *) d;
 397
 398     /* Starting point */
 399     x = gsl_vector_alloc (my_func.n);
 400     for (i = 0; (i < my_func.n); i++)
 401     {
 402         gsl_vector_set (x, i, d->params[i]);
 403     }
 404
 405     /* Step size, different for each of the parameters */
 406     dx = gsl_vector_alloc (my_func.n);
 407     for (i = 0; (i < my_func.n); i++)
 408     {
 409         gsl_vector_set (dx, i, 0.1*d->params[i]);
 410     }
 411
 412     T = gsl_multimin_fminimizer_nmsimplex;
 413     s = gsl_multimin_fminimizer_alloc (T, my_func.n);
 414
 415     gsl_multimin_fminimizer_set (s, &my_func, x, dx);
 416     gsl_vector_free (x);
 417     gsl_vector_free (dx);
 418
 419     printf ("%5s", "Iter");
 420     for (i = 0; (i < my_func.n); i++)
 421     {
 422         printf(" %12s", epnm(my_func.n, i));
 423     }
 424     printf (" %12s %12s\n", "NM Size", "Chi2");
 425
 426     do
 427     {
 428         iter++;
 429         status = gsl_multimin_fminimizer_iterate (s);
 430
 431         if (status != 0)
 432         {
 433             gmx_fatal(FARGS, "Something went wrong in the iteration in minimizer %s",
 434                       gsl_multimin_fminimizer_name(s));
 435         }
 436
 437         d2     = gsl_multimin_fminimizer_minimum(s);
 438         size   = gsl_multimin_fminimizer_size(s);
 439         status = gsl_multimin_test_size(size, tol);
 440
 441         if (status == GSL_SUCCESS)
 442         {
 443             printf ("Minimum found using %s at:\n",
 444                     gsl_multimin_fminimizer_name(s));
 445         }
 446
 447         printf ("%5d", iter);
 448         for (i = 0; (i < my_func.n); i++)
 449         {
 450             printf(" %12.4e", gsl_vector_get (s->x, i));
 451         }
 452         printf (" %12.4e %12.4e\n", size, d2);
 453     }
 454     while ((status == GSL_CONTINUE) && (iter < maxiter));
 455
 456     gsl_multimin_fminimizer_free (s);
 457 }
 458
 459 static void preprocess_remd(FILE *fp, t_remd_data *d, real cutoff, real tref,
 460                             real ucut, gmx_bool bBack, real Euf, real Efu,
 461                             real Ei, real t0, real t1, gmx_bool bSum, gmx_bool bDiscrete,
 462                             int nmult)
 463 {
 464     int  i, j, ninter;
 465     real dd, tau_f, tau_u;
 466
 467     ninter = (ucut > cutoff) ? 1 : 0;
 468     if (ninter && (ucut <= cutoff))
 469     {
 470         gmx_fatal(FARGS, "You have requested an intermediate but the cutoff for intermediates %f is smaller than the normal cutoff(%f)", ucut, cutoff);
 471     }
 472
 473     if (!bBack)
 474     {
 475         d->nparams = 2;
 476         d->nstate  = 2;
 477     }
 478     else
 479     {
 480         d->nparams = 4*(1+ninter);
 481         d->nstate  = 2+ninter;
 482     }
 483     d->bSum      = bSum;
 484     d->bDiscrete = bDiscrete;
 485     snew(d->beta, d->nreplica);
 486     snew(d->state, d->nreplica);
 487     snew(d->bMask, d->nreplica);
 488     snew(d->d2_replica, d->nreplica);
 489     snew(d->sumft, d->nframe);
 490     snew(d->sumit, d->nframe);
 491     snew(d->sumfct, d->nframe);
 492     snew(d->sumict, d->nframe);
 493     snew(d->params, d->nparams);
 494     snew(d->fcalt, d->nreplica);
 495     snew(d->icalt, d->nreplica);
 496
 497     /* convert_times(d->nframe,d->time); */
 498
 499     if (t0 < 0)
 500     {
 501         d->j0 = 0;
 502     }
 503     else
 504     {
 505         for (d->j0 = 0; (d->j0 < d->nframe) && (d->time[d->j0] < t0); d->j0++)
 506         {
 507             ;
 508         }
 509     }
 510     if (t1 < 0)
 511     {
 512         d->j1 = d->nframe;
 513     }
 514     else
 515     {
 516         for (d->j1 = 0; (d->j1 < d->nframe) && (d->time[d->j1] < t1); d->j1++)
 517         {
 518             ;
 519         }
 520     }
 521     if ((d->j1-d->j0) < d->nparams+2)
 522     {
 523         gmx_fatal(FARGS, "Start (%f) or end time (%f) for fitting inconsistent. Reduce t0, increase t1 or supply more data", t0, t1);
 524     }
 525     fprintf(fp, "Will optimize from %g to %g\n",
 526             d->time[d->j0], d->time[d->j1-1]);
 527     d->nmask = d->nreplica;
 528     for (i = 0; (i < d->nreplica); i++)
 529     {
 530         snew(d->beta[i], d->nframe);
 531         snew(d->state[i], d->nframe);
 532         snew(d->fcalt[i], d->nframe);
 533         snew(d->icalt[i], d->nframe);
 534         d->bMask[i] = TRUE;
 535         for (j = 0; (j < d->nframe); j++)
 536         {
 537             d->beta[i][j] = 1.0/(BOLTZ*d->temp[i][j]);
 538             dd            = d->data[i][j];
 539             if (bDiscrete)
 540             {
 541                 if (dd <= cutoff)
 542                 {
 543                     d->state[i][j] = 0;
 544                 }
 545                 else if ((ucut > cutoff) && (dd <= ucut))
 546                 {
 547                     d->state[i][j] = 1;
 548                 }
 549                 else
 550                 {
 551                     d->state[i][j] = d->nstate-1;
 552                 }
 553             }
 554             else
 555             {
 556                 d->state[i][j] = dd*nmult;
 557             }
 558         }
 559     }
 560     sum_ft(d);
 561
 562     /* Assume forward rate constant is half the total time in this
 563      * simulation and backward is ten times as long */
 564     if (bDiscrete)
 565     {
 566         tau_f            = d->time[d->nframe-1];
 567         tau_u            = 4*tau_f;
 568         d->params[epEuf] = Euf;
 569         d->params[epAuf] = exp(d->params[epEuf]/(BOLTZ*tref))/tau_f;
 570         if (bBack)
 571         {
 572             d->params[epEfu] = Efu;
 573             d->params[epAfu] = exp(d->params[epEfu]/(BOLTZ*tref))/tau_u;
 574             if (ninter > 0)
 575             {
 576                 d->params[eqEui] = Ei;
 577                 d->params[eqAui] = exp(d->params[eqEui]/(BOLTZ*tref))/tau_u;
 578                 d->params[eqEiu] = Ei;
 579                 d->params[eqAiu] = exp(d->params[eqEiu]/(BOLTZ*tref))/tau_u;
 580             }
 581         }
 582         else
 583         {
 584             d->params[epAfu]  = 0;
 585             d->params[epEfu]  = 0;
 586         }
 587     }
 588     else
 589     {
 590         d->params[epEuf] = Euf;
 591         if (d->data2)
 592         {
 593             d->params[epAuf] = 0.1;
 594         }
 595         else
 596         {
 597             d->params[epAuf] = 20.0;
 598         }
 599     }
 600 }
 601
 602 static real tau(real A, real E, real T)
 603 {
 604     return exp(E/(BOLTZ*T))/A;
 605 }
 606
 607 static real folded_fraction(t_remd_data *d, real tref)
 608 {
 609     real tauf, taub;
 610
 611     tauf = tau(d->params[epAuf], d->params[epEuf], tref);
 612     taub = tau(d->params[epAfu], d->params[epEfu], tref);
 613
 614     return (taub/(tauf+taub));
 615 }
 616
 617 static void print_tau(FILE *gp, t_remd_data *d, real tref)
 618 {
 619     real tauf, taub, ddd, fff, DG, DH, TDS, Tm, Tb, Te, Fb, Fe, Fm;
 620     int  i, np = d->nparams;
 621
 622     ddd = calc_d2(d);
 623     fprintf(gp, "Final value for Chi2 = %12.5e (%d replicas)\n", ddd, d->nmask);
 624     tauf = tau(d->params[epAuf], d->params[epEuf], tref);
 625     fprintf(gp, "%s = %12.5e %s = %12.5e (kJ/mole)\n",
 626             epnm(np, epAuf), d->params[epAuf],
 627             epnm(np, epEuf), d->params[epEuf]);
 628     if (bBack(d))
 629     {
 630         taub = tau(d->params[epAfu], d->params[epEfu], tref);
 631         fprintf(gp, "%s = %12.5e %s = %12.5e (kJ/mole)\n",
 632                 epnm(np, epAfu), d->params[epAfu],
 633                 epnm(np, epEfu), d->params[epEfu]);
 634         fprintf(gp, "Equilibrium properties at T = %g\n", tref);
 635         fprintf(gp, "tau_f = %8.3f ns, tau_b = %8.3f ns\n", tauf/1000, taub/1000);
 636         fff = taub/(tauf+taub);
 637         DG  = BOLTZ*tref*log(fff/(1-fff));
 638         DH  = d->params[epEfu]-d->params[epEuf];
 639         TDS = DH-DG;
 640         fprintf(gp, "Folded fraction     F = %8.3f\n", fff);
 641         fprintf(gp, "Unfolding energies: DG = %8.3f  DH = %8.3f TDS = %8.3f\n",
 642                 DG, DH, TDS);
 643         Tb = 260;
 644         Te = 420;
 645         Tm = 0;
 646         Fm = 0;
 647         Fb = folded_fraction(d, Tb);
 648         Fe = folded_fraction(d, Te);
 649         while ((Te-Tb > 0.001) && (Fm != 0.5))
 650         {
 651             Tm = 0.5*(Tb+Te);
 652             Fm = folded_fraction(d, Tm);
 653             if (Fm > 0.5)
 654             {
 655                 Fb = Fm;
 656                 Tb = Tm;
 657             }
 658             else if (Fm < 0.5)
 659             {
 660                 Te = Tm;
 661                 Fe = Fm;
 662             }
 663         }
 664         if ((Fb-0.5)*(Fe-0.5) <= 0)
 665         {
 666             fprintf(gp, "Melting temperature Tm = %8.3f K\n", Tm);
 667         }
 668         else
 669         {
 670             fprintf(gp, "No melting temperature detected between 260 and 420K\n");
 671         }
 672         if (np > 4)
 673         {
 674             char *ptr;
 675             fprintf(gp, "Data for intermediates at T = %g\n", tref);
 676             fprintf(gp, "%8s  %10s  %10s  %10s\n", "Name", "A", "E", "tau");
 677             for (i = 0; (i < np/2); i++)
 678             {
 679                 tauf = tau(d->params[2*i], d->params[2*i+1], tref);
 680                 ptr  = epnm(d->nparams, 2*i);
 681                 fprintf(gp, "%8s  %10.3e  %10.3e  %10.3e\n", ptr+1,
 682                         d->params[2*i], d->params[2*i+1], tauf/1000);
 683             }
 684         }
 685     }
 686     else
 687     {
 688         fprintf(gp, "Equilibrium properties at T = %g\n", tref);
 689         fprintf(gp, "tau_f = %8.3f\n", tauf);
 690     }
 691 }
 692
 693 static void dump_remd_parameters(FILE *gp, t_remd_data *d, const char *fn,
 694                                  const char *fn2, const char *rfn,
 695                                  const char *efn, const char *mfn, int skip, real tref,
 696                                  output_env_t oenv)
 697 {
 698     FILE       *fp, *hp;
 699     int         i, j, np = d->nparams;
 700     real        rhs, tauf, taub, fff, DG;
 701     real       *params;
 702     const char *leg[]  = { "Measured", "Fit", "Difference" };
 703     const char *mleg[] = { "Folded fraction", "DG (kJ/mole)"};
 704     char      **rleg;
 705     real        fac[] = { 0.97, 0.98, 0.99, 1.0, 1.01, 1.02, 1.03 };
 706 #define NFAC asize(fac)
 707     real        d2[NFAC];
 708     double      norm;
 709
 710     integrate_dfdt(d);
 711     print_tau(gp, d, tref);
 712     norm = (d->bDiscrete ? 1.0/d->nmask : 1.0);
 713
 714     if (fn)
 715     {
 716         fp = xvgropen(fn, "Optimized fit to data", "Time (ps)", "Fraction Folded", oenv);
 717         xvgr_legend(fp, asize(leg), leg, oenv);
 718         for (i = 0; (i < d->nframe); i++)
 719         {
 720             if ((skip <= 0) || ((i % skip) == 0))
 721             {
 722                 fprintf(fp, "%12.5e  %12.5e  %12.5e  %12.5e\n", d->time[i],
 723                         d->sumft[i]*norm, d->sumfct[i]*norm,
 724                         (d->sumft[i]-d->sumfct[i])*norm);
 725             }
 726         }
 727         ffclose(fp);
 728     }
 729     if (!d->bSum && rfn)
 730     {
 731         snew(rleg, d->nreplica*2);
 732         for (i = 0; (i < d->nreplica); i++)
 733         {
 734             snew(rleg[2*i], 32);
 735             snew(rleg[2*i+1], 32);
 736             sprintf(rleg[2*i], "\\f{4}F(t) %d", i);
 737             sprintf(rleg[2*i+1], "\\f{12}F \\f{4}(t) %d", i);
 738         }
 739         fp = xvgropen(rfn, "Optimized fit to data", "Time (ps)", "Fraction Folded", oenv);
 740         xvgr_legend(fp, d->nreplica*2, (const char**)rleg, oenv);
 741         for (j = 0; (j < d->nframe); j++)
 742         {
 743             if ((skip <= 0) || ((j % skip) == 0))
 744             {
 745                 fprintf(fp, "%12.5e", d->time[j]);
 746                 for (i = 0; (i < d->nreplica); i++)
 747                 {
 748                     fprintf(fp, "  %5f  %9.2e", is_folded(d, i, j), d->fcalt[i][j]);
 749                 }
 750                 fprintf(fp, "\n");
 751             }
 752         }
 753         ffclose(fp);
 754     }
 755
 756     if (fn2 && (d->nstate > 2))
 757     {
 758         fp = xvgropen(fn2, "Optimized fit to data", "Time (ps)",
 759                       "Fraction Intermediate", oenv);
 760         xvgr_legend(fp, asize(leg), leg, oenv);
 761         for (i = 0; (i < d->nframe); i++)
 762         {
 763             if ((skip <= 0) || ((i % skip) == 0))
 764             {
 765                 fprintf(fp, "%12.5e  %12.5e  %12.5e  %12.5e\n", d->time[i],
 766                         d->sumit[i]*norm, d->sumict[i]*norm,
 767                         (d->sumit[i]-d->sumict[i])*norm);
 768             }
 769         }
 770         ffclose(fp);
 771     }
 772     if (mfn)
 773     {
 774         if (bBack(d))
 775         {
 776             fp = xvgropen(mfn, "Melting curve", "T (K)", "", oenv);
 777             xvgr_legend(fp, asize(mleg), mleg, oenv);
 778             for (i = 260; (i <= 420); i++)
 779             {
 780                 tauf = tau(d->params[epAuf], d->params[epEuf], 1.0*i);
 781                 taub = tau(d->params[epAfu], d->params[epEfu], 1.0*i);
 782                 fff  = taub/(tauf+taub);
 783                 DG   = BOLTZ*i*log(fff/(1-fff));
 784                 fprintf(fp, "%5d  %8.3f  %8.3f\n", i, fff, DG);
 785             }
 786             ffclose(fp);
 787         }
 788     }
 789
 790     if (efn)
 791     {
 792         snew(params, d->nparams);
 793         for (i = 0; (i < d->nparams); i++)
 794         {
 795             params[i] = d->params[i];
 796         }
 797
 798         hp = xvgropen(efn, "Chi2 as a function of relative parameter",
 799                       "Fraction", "Chi2", oenv);
 800         for (j = 0; (j < d->nparams); j++)
 801         {
 802             /* Reset all parameters to optimized values */
 803             fprintf(hp, "@type xy\n");
 804             for (i = 0; (i < d->nparams); i++)
 805             {
 806                 d->params[i] = params[i];
 807             }
 808             /* Now modify one of them */
 809             for (i = 0; (i < NFAC); i++)
 810             {
 811                 d->params[j] = fac[i]*params[j];
 812                 d2[i]        = calc_d2(d);
 813                 fprintf(gp, "%s = %12g  d2 = %12g\n", epnm(np, j), d->params[j], d2[i]);
 814                 fprintf(hp, "%12g  %12g\n", fac[i], d2[i]);
 815             }
 816             fprintf(hp, "&\n");
 817         }
 818         ffclose(hp);
 819         for (i = 0; (i < d->nparams); i++)
 820         {
 821             d->params[i] = params[i];
 822         }
 823         sfree(params);
 824     }
 825     if (!d->bSum)
 826     {
 827         for (i = 0; (i < d->nreplica); i++)
 828         {
 829             fprintf(gp, "Chi2[%3d] = %8.2e\n", i, d->d2_replica[i]);
 830         }
 831     }
 832 }
 833 #endif /*HAVE_LIBGSL*/
 834
 835 int gmx_kinetics(int argc, char *argv[])
 836 {
 837     const char     *desc[] = {
 838         "[TT]g_kinetics[tt] reads two [TT].xvg[tt] files, each one containing data for N replicas.",
 839         "The first file contains the temperature of each replica at each timestep,",
 840         "and the second contains real values that can be interpreted as",
 841         "an indicator for folding. If the value in the file is larger than",
 842         "the cutoff it is taken to be unfolded and the other way around.[PAR]",
 843         "From these data an estimate of the forward and backward rate constants",
 844         "for folding is made at a reference temperature. In addition,",
 845         "a theoretical melting curve and free energy as a function of temperature",
 846         "are printed in an [TT].xvg[tt] file.[PAR]",
 847         "The user can give a max value to be regarded as intermediate",
 848         "([TT]-ucut[tt]), which, when given will trigger the use of an intermediate state",
 849         "in the algorithm to be defined as those structures that have",
 850         "cutoff < DATA < ucut. Structures with DATA values larger than ucut will",
 851         "not be regarded as potential folders. In this case 8 parameters are optimized.[PAR]",
 852         "The average fraction foled is printed in an [TT].xvg[tt] file together with the fit to it.",
 853         "If an intermediate is used a further file will show the build of the intermediate and the fit to that process.[PAR]",
 854         "The program can also be used with continuous variables (by setting",
 855         "[TT]-nodiscrete[tt]). In this case kinetics of other processes can be",
 856         "studied. This is very much a work in progress and hence the manual",
 857         "(this information) is lagging behind somewhat.[PAR]",
 858         "In order to compile this program you need access to the GNU",
 859         "scientific library."
 860     };
 861     static int      nreplica  = 1;
 862     static real     tref      = 298.15;
 863     static real     cutoff    = 0.2;
 864     static real     ucut      = 0.0;
 865     static real     Euf       = 10;
 866     static real     Efu       = 30;
 867     static real     Ei        = 10;
 868     static gmx_bool bHaveT    = TRUE;
 869     static real     t0        = -1;
 870     static real     t1        = -1;
 871     static real     tb        = 0;
 872     static real     te        = 0;
 873     static real     tol       = 1e-3;
 874     static int      maxiter   = 100;
 875     static int      skip      = 0;
 876     static int      nmult     = 1;
 877     static gmx_bool bBack     = TRUE;
 878     static gmx_bool bSplit    = TRUE;
 879     static gmx_bool bSum      = TRUE;
 880     static gmx_bool bDiscrete = TRUE;
 881     t_pargs         pa[]      = {
 882         { "-time",    FALSE, etBOOL, {&bHaveT},
 883           "Expect a time in the input" },
 884         { "-b",       FALSE, etREAL, {&tb},
 885           "First time to read from set" },
 886         { "-e",       FALSE, etREAL, {&te},
 887           "Last time to read from set" },
 888         { "-bfit",    FALSE, etREAL, {&t0},
 889           "Time to start the fit from" },
 890         { "-efit",    FALSE, etREAL, {&t1},
 891           "Time to end the fit" },
 892         { "-T",       FALSE, etREAL, {&tref},
 893           "Reference temperature for computing rate constants" },
 894         { "-n",       FALSE, etINT, {&nreplica},
 895           "Read data for this number of replicas. Only necessary when files are written in xmgrace format using @type and & as delimiters." },
 896         { "-cut",     FALSE, etREAL, {&cutoff},
 897           "Cut-off (max) value for regarding a structure as folded" },
 898         { "-ucut",    FALSE, etREAL, {&ucut},
 899           "Cut-off (max) value for regarding a structure as intermediate (if not folded)" },
 900         { "-euf",     FALSE, etREAL, {&Euf},
 901           "Initial guess for energy of activation for folding (kJ/mol)" },
 902         { "-efu",     FALSE, etREAL, {&Efu},
 903           "Initial guess for energy of activation for unfolding (kJ/mol)" },
 904         { "-ei",      FALSE, etREAL, {&Ei},
 905           "Initial guess for energy of activation for intermediates (kJ/mol)" },
 906         { "-maxiter", FALSE, etINT, {&maxiter},
 907           "Max number of iterations" },
 908         { "-back",    FALSE, etBOOL, {&bBack},
 909           "Take the back reaction into account" },
 910         { "-tol",     FALSE, etREAL, {&tol},
 911           "Absolute tolerance for convergence of the Nelder and Mead simplex algorithm" },
 912         { "-skip",    FALSE, etINT, {&skip},
 913           "Skip points in the output [TT].xvg[tt] file" },
 914         { "-split",   FALSE, etBOOL, {&bSplit},
 915           "Estimate error by splitting the number of replicas in two and refitting" },
 916         { "-sum",     FALSE, etBOOL, {&bSum},
 917           "Average folding before computing [GRK]chi[grk]^2" },
 918         { "-discrete", FALSE, etBOOL, {&bDiscrete},
 919           "Use a discrete folding criterion (F <-> U) or a continuous one" },
 920         { "-mult",    FALSE, etINT, {&nmult},
 921           "Factor to multiply the data with before discretization" }
 922     };
 923 #define NPA asize(pa)
 924
 925     FILE        *fp;
 926     real         dt_t, dt_d, dt_d2;
 927     int          nset_t, nset_d, nset_d2, n_t, n_d, n_d2, i;
 928     const char  *tfile, *dfile, *dfile2;
 929     t_remd_data  remd;
 930     output_env_t oenv;
 931
 932     t_filenm     fnm[] = {
 933         { efXVG, "-f",    "temp",    ffREAD   },
 934         { efXVG, "-d",    "data",    ffREAD   },
 935         { efXVG, "-d2",   "data2",   ffOPTRD  },
 936         { efXVG, "-o",    "ft_all",  ffWRITE  },
 937         { efXVG, "-o2",   "it_all",  ffOPTWR  },
 938         { efXVG, "-o3",   "ft_repl", ffOPTWR  },
 939         { efXVG, "-ee",   "err_est", ffOPTWR  },
 940         { efLOG, "-g",    "remd",    ffWRITE  },
 941         { efXVG, "-m",    "melt",    ffWRITE  }
 942     };
 943 #define NFILE asize(fnm)
 944
 945     parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_BE_NICE | PCA_TIME_UNIT,
 946                       NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv);
 947
 948 #ifdef HAVE_LIBGSL
 949     please_cite(stdout, "Spoel2006d");
 950     if (cutoff < 0)
 951     {
 952         gmx_fatal(FARGS, "cutoff should be >= 0 (rather than %f)", cutoff);
 953     }
 954
 955     tfile   = opt2fn("-f", NFILE, fnm);
 956     dfile   = opt2fn("-d", NFILE, fnm);
 957     dfile2  = opt2fn_null("-d2", NFILE, fnm);
 958
 959     fp = ffopen(opt2fn("-g", NFILE, fnm), "w");
 960
 961     remd.temp = read_xvg_time(tfile, bHaveT,
 962                               opt2parg_bSet("-b", NPA, pa), tb,
 963                               opt2parg_bSet("-e", NPA, pa), te,
 964                               nreplica, &nset_t, &n_t, &dt_t, &remd.time);
 965     printf("Read %d sets of %d points in %s, dt = %g\n\n", nset_t, n_t, tfile, dt_t);
 966     sfree(remd.time);
 967
 968     remd.data = read_xvg_time(dfile, bHaveT,
 969                               opt2parg_bSet("-b", NPA, pa), tb,
 970                               opt2parg_bSet("-e", NPA, pa), te,
 971                               nreplica, &nset_d, &n_d, &dt_d, &remd.time);
 972     printf("Read %d sets of %d points in %s, dt = %g\n\n", nset_d, n_d, dfile, dt_d);
 973
 974     if ((nset_t != nset_d) || (n_t != n_d) || (dt_t != dt_d))
 975     {
 976         gmx_fatal(FARGS, "Files %s and %s are inconsistent. Check log file",
 977                   tfile, dfile);
 978     }
 979
 980     if (dfile2)
 981     {
 982         remd.data2 = read_xvg_time(dfile2, bHaveT,
 983                                    opt2parg_bSet("-b", NPA, pa), tb,
 984                                    opt2parg_bSet("-e", NPA, pa), te,
 985                                    nreplica, &nset_d2, &n_d2, &dt_d2, &remd.time);
 986         printf("Read %d sets of %d points in %s, dt = %g\n\n",
 987                nset_d2, n_d2, dfile2, dt_d2);
 988         if ((nset_d2 != nset_d) || (n_d != n_d2) || (dt_d != dt_d2))
 989         {
 990             gmx_fatal(FARGS, "Files %s and %s are inconsistent. Check log file",
 991                       dfile, dfile2);
 992         }
 993     }
 994     else
 995     {
 996         remd.data2 = NULL;
 997     }
 998
 999     remd.nreplica  = nset_d;
1000     remd.nframe    = n_d;
1001     remd.dt        = 1;
1002     preprocess_remd(fp, &remd, cutoff, tref, ucut, bBack, Euf, Efu, Ei, t0, t1,
1003                     bSum, bDiscrete, nmult);
1004
1005     optimize_remd_parameters(&remd, maxiter, tol);
1006
1007     dump_remd_parameters(fp, &remd, opt2fn("-o", NFILE, fnm),
1008                          opt2fn_null("-o2", NFILE, fnm),
1009                          opt2fn_null("-o3", NFILE, fnm),
1010                          opt2fn_null("-ee", NFILE, fnm),
1011                          opt2fn("-m", NFILE, fnm), skip, tref, oenv);
1012
1013     if (bSplit)
1014     {
1015         printf("Splitting set of replicas in two halves\n");
1016         for (i = 0; (i < remd.nreplica); i++)
1017         {
1018             remd.bMask[i] = FALSE;
1019         }
1020         remd.nmask = 0;
1021         for (i = 0; (i < remd.nreplica); i += 2)
1022         {
1023             remd.bMask[i] = TRUE;
1024             remd.nmask++;
1025         }
1026         sum_ft(&remd);
1027         optimize_remd_parameters(&remd, maxiter, tol);
1028         dump_remd_parameters(fp, &remd, "test1.xvg", NULL, NULL, NULL, NULL, skip, tref, oenv);
1029
1030         for (i = 0; (i < remd.nreplica); i++)
1031         {
1032             remd.bMask[i] = !remd.bMask[i];
1033         }
1034         remd.nmask = remd.nreplica - remd.nmask;
1035
1036         sum_ft(&remd);
1037         optimize_remd_parameters(&remd, maxiter, tol);
1038         dump_remd_parameters(fp, &remd, "test2.xvg", NULL, NULL, NULL, NULL, skip, tref, oenv);
1039
1040         for (i = 0; (i < remd.nreplica); i++)
1041         {
1042             remd.bMask[i] = FALSE;
1043         }
1044         remd.nmask = 0;
1045         for (i = 0; (i < remd.nreplica/2); i++)
1046         {
1047             remd.bMask[i] = TRUE;
1048             remd.nmask++;
1049         }
1050         sum_ft(&remd);
1051         optimize_remd_parameters(&remd, maxiter, tol);
1052         dump_remd_parameters(fp, &remd, "test1.xvg", NULL, NULL, NULL, NULL, skip, tref, oenv);
1053
1054         for (i = 0; (i < remd.nreplica); i++)
1055         {
1056             remd.bMask[i] = FALSE;
1057         }
1058         remd.nmask = 0;
1059         for (i = remd.nreplica/2; (i < remd.nreplica); i++)
1060         {
1061             remd.bMask[i] = TRUE;
1062             remd.nmask++;
1063         }
1064         sum_ft(&remd);
1065         optimize_remd_parameters(&remd, maxiter, tol);
1066         dump_remd_parameters(fp, &remd, "test1.xvg", NULL, NULL, NULL, NULL, skip, tref, oenv);
1067     }
1068     ffclose(fp);
1069
1070     view_all(oenv, NFILE, fnm);
1071
1072 #else
1073     fprintf(stderr, "This program should be compiled with the GNU scientific library. Please install the library and reinstall GROMACS.\n");
1074 #endif /*HAVE_LIBGSL*/
1075
1076     return 0;
1077 }