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