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