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