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 #ifdef HAVE_LIBGSL
  59 #include <gsl/gsl_multimin.h>
  60
  61 enum { epAuf, epEuf, epAfu, epEfu, epNR };
  62 enum { eqAif, eqEif, eqAfi, eqEfi, eqAui, eqEui, eqAiu, eqEiu, eqNR };
  63 static char *eep[epNR] = { "Af", "Ef", "Au", "Eu" };
  64 static char *eeq[eqNR] = { "Aif","Eif","Afi","Efi","Aui","Eui","Aiu","Eiu" };
  65
  66 typedef struct {
  67   int nreplica;   /* Number of replicas in the calculation                   */
  68   int nframe;     /* Number of time frames                                   */
  69   int nstate;     /* Number of states the system can be in, e.g. F,I,U       */
  70   int nparams;    /* Is 2, 4 or 8                                            */
  71   gmx_bool *bMask;    /* Determine whether this replica is part of the d2 comp.  */
  72   gmx_bool bSum;
  73   gmx_bool bDiscrete; /* Use either discrete folding (0/1) or a continuous       */
  74                   /* criterion */
  75   int nmask;      /* Number of replicas taken into account                   */
  76   real dt;        /* Timestep between frames                                 */
  77   int  j0,j1;     /* Range of frames used in calculating delta               */
  78   real **temp,**data,**data2;
  79   int  **state;   /* State index running from 0 (F) to nstate-1 (U)          */
  80   real **beta,**fcalt,**icalt;
  81   real *time,*sumft,*sumit,*sumfct,*sumict;
  82   real *params;
  83   real *d2_replica;
  84 } t_remd_data;
  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
 668 int gmx_kinetics(int argc,char *argv[])
 669 {
 670   const char *desc[] = {
 671     "[TT]g_kinetics[tt] reads two [TT].xvg[tt] files, each one containing data for N replicas.",
 672     "The first file contains the temperature of each replica at each timestep,",
 673     "and the second contains real values that can be interpreted as",
 674     "an indicator for folding. If the value in the file is larger than",
 675     "the cutoff it is taken to be unfolded and the other way around.[PAR]",
 676     "From these data an estimate of the forward and backward rate constants",
 677     "for folding is made at a reference temperature. In addition,",
 678     "a theoretical melting curve and free energy as a function of temperature",
 679     "are printed in an [TT].xvg[tt] file.[PAR]",
 680     "The user can give a max value to be regarded as intermediate",
 681     "([TT]-ucut[tt]), which, when given will trigger the use of an intermediate state",
 682     "in the algorithm to be defined as those structures that have",
 683     "cutoff < DATA < ucut. Structures with DATA values larger than ucut will",
 684     "not be regarded as potential folders. In this case 8 parameters are optimized.[PAR]",
 685     "The average fraction foled is printed in an [TT].xvg[tt] file together with the fit to it.",
 686     "If an intermediate is used a further file will show the build of the intermediate and the fit to that process.[PAR]",
 687     "The program can also be used with continuous variables (by setting",
 688     "[TT]-nodiscrete[tt]). In this case kinetics of other processes can be",
 689     "studied. This is very much a work in progress and hence the manual",
 690     "(this information) is lagging behind somewhat.[PAR]",
 691     "In order to compile this program you need access to the GNU",
 692     "scientific library."
 693   };
 694   static int  nreplica  = 1;
 695   static real tref      = 298.15;
 696   static real cutoff    = 0.2;
 697   static real ucut      = 0.0;
 698   static real Euf       = 10;
 699   static real Efu       = 30;
 700   static real Ei        = 10;
 701   static gmx_bool bHaveT    = TRUE;
 702   static real t0        = -1;
 703   static real t1        = -1;
 704   static real tb        = 0;
 705   static real te        = 0;
 706   static real tol       = 1e-3;
 707   static int  maxiter   = 100;
 708   static int  skip      = 0;
 709   static int  nmult     = 1;
 710   static gmx_bool bBack     = TRUE;
 711   static gmx_bool bSplit    = TRUE;
 712   static gmx_bool bSum      = TRUE;
 713   static gmx_bool bDiscrete = TRUE;
 714   t_pargs pa[] = {
 715     { "-time",    FALSE, etBOOL, {&bHaveT},
 716       "Expect a time in the input" },
 717     { "-b",       FALSE, etREAL, {&tb},
 718       "First time to read from set" },
 719     { "-e",       FALSE, etREAL, {&te},
 720       "Last time to read from set" },
 721     { "-bfit",    FALSE, etREAL, {&t0},
 722       "Time to start the fit from" },
 723     { "-efit",    FALSE, etREAL, {&t1},
 724       "Time to end the fit" },
 725     { "-T",       FALSE, etREAL, {&tref},
 726       "Reference temperature for computing rate constants" },
 727     { "-n",       FALSE, etINT, {&nreplica},
 728       "Read data for # replicas. Only necessary when files are written in xmgrace format using @type and & as delimiters." },
 729     { "-cut",     FALSE, etREAL, {&cutoff},
 730       "Cut-off (max) value for regarding a structure as folded" },
 731     { "-ucut",    FALSE, etREAL, {&ucut},
 732       "Cut-off (max) value for regarding a structure as intermediate (if not folded)" },
 733     { "-euf",     FALSE, etREAL, {&Euf},
 734       "Initial guess for energy of activation for folding (kJ/mole)" },
 735     { "-efu",     FALSE, etREAL, {&Efu},
 736       "Initial guess for energy of activation for unfolding (kJ/mole)" },
 737     { "-ei",      FALSE, etREAL, {&Ei},
 738       "Initial guess for energy of activation for intermediates (kJ/mole)" },
 739     { "-maxiter", FALSE, etINT, {&maxiter},
 740       "Max number of iterations" },
 741     { "-back",    FALSE, etBOOL, {&bBack},
 742       "Take the back reaction into account" },
 743     { "-tol",     FALSE, etREAL,{&tol},
 744       "Absolute tolerance for convergence of the Nelder and Mead simplex algorithm" },
 745     { "-skip",    FALSE, etINT, {&skip},
 746       "Skip points in the output [TT].xvg[tt] file" },
 747     { "-split",   FALSE, etBOOL,{&bSplit},
 748       "Estimate error by splitting the number of replicas in two and refitting" },
 749     { "-sum",     FALSE, etBOOL,{&bSum},
 750       "Average folding before computing chi^2" },
 751     { "-discrete", FALSE, etBOOL, {&bDiscrete},
 752       "Use a discrete folding criterium (F <-> U) or a continuous one" },
 753     { "-mult",    FALSE, etINT, {&nmult},
 754       "Factor to multiply the data with before discretization" }
 755   };
 756 #define NPA asize(pa)
 757
 758   FILE        *fp;
 759   real        dt_t,dt_d,dt_d2;
 760   int         nset_t,nset_d,nset_d2,n_t,n_d,n_d2,i;
 761   const char  *tfile,*dfile,*dfile2;
 762   t_remd_data remd;
 763   output_env_t oenv;
 764
 765   t_filenm fnm[] = {
 766     { efXVG, "-f",    "temp",    ffREAD   },
 767     { efXVG, "-d",    "data",    ffREAD   },
 768     { efXVG, "-d2",   "data2",   ffOPTRD  },
 769     { efXVG, "-o",    "ft_all",  ffWRITE  },
 770     { efXVG, "-o2",   "it_all",  ffOPTWR  },
 771     { efXVG, "-o3",   "ft_repl", ffOPTWR  },
 772     { efXVG, "-ee",   "err_est", ffOPTWR  },
 773     { efLOG, "-g",    "remd",    ffWRITE  },
 774     { efXVG, "-m",    "melt",    ffWRITE  }
 775   };
 776 #define NFILE asize(fnm)
 777
 778   CopyRight(stderr,argv[0]);
 779   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_BE_NICE | PCA_TIME_UNIT,
 780                     NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL,&oenv);
 781
 782   please_cite(stdout,"Spoel2006d");
 783   if (cutoff < 0)
 784     gmx_fatal(FARGS,"cutoff should be >= 0 (rather than %f)",cutoff);
 785
 786   tfile   = opt2fn("-f",NFILE,fnm);
 787   dfile   = opt2fn("-d",NFILE,fnm);
 788   dfile2  = opt2fn_null("-d2",NFILE,fnm);
 789
 790   fp = ffopen(opt2fn("-g",NFILE,fnm),"w");
 791
 792   remd.temp = read_xvg_time(tfile,bHaveT,
 793                             opt2parg_bSet("-b",NPA,pa),tb,
 794                             opt2parg_bSet("-e",NPA,pa),te,
 795                             nreplica,&nset_t,&n_t,&dt_t,&remd.time);
 796   printf("Read %d sets of %d points in %s, dt = %g\n\n",nset_t,n_t,tfile,dt_t);
 797   sfree(remd.time);
 798
 799   remd.data = read_xvg_time(dfile,bHaveT,
 800                             opt2parg_bSet("-b",NPA,pa),tb,
 801                             opt2parg_bSet("-e",NPA,pa),te,
 802                             nreplica,&nset_d,&n_d,&dt_d,&remd.time);
 803   printf("Read %d sets of %d points in %s, dt = %g\n\n",nset_d,n_d,dfile,dt_d);
 804
 805   if ((nset_t != nset_d) || (n_t != n_d) || (dt_t != dt_d))
 806     gmx_fatal(FARGS,"Files %s and %s are inconsistent. Check log file",
 807               tfile,dfile);
 808
 809   if (dfile2) {
 810     remd.data2 = read_xvg_time(dfile2,bHaveT,
 811                                opt2parg_bSet("-b",NPA,pa),tb,
 812                                opt2parg_bSet("-e",NPA,pa),te,
 813                                nreplica,&nset_d2,&n_d2,&dt_d2,&remd.time);
 814     printf("Read %d sets of %d points in %s, dt = %g\n\n",
 815            nset_d2,n_d2,dfile2,dt_d2);
 816     if ((nset_d2 != nset_d) || (n_d != n_d2) || (dt_d != dt_d2))
 817       gmx_fatal(FARGS,"Files %s and %s are inconsistent. Check log file",
 818                 dfile,dfile2);
 819   }
 820   else {
 821     remd.data2 = NULL;
 822   }
 823
 824   remd.nreplica  = nset_d;
 825   remd.nframe    = n_d;
 826   remd.dt        = 1;
 827   preprocess_remd(fp,&remd,cutoff,tref,ucut,bBack,Euf,Efu,Ei,t0,t1,
 828                   bSum,bDiscrete,nmult);
 829
 830   optimize_remd_parameters(fp,&remd,maxiter,tol);
 831
 832   dump_remd_parameters(fp,&remd,opt2fn("-o",NFILE,fnm),
 833                        opt2fn_null("-o2",NFILE,fnm),
 834                        opt2fn_null("-o3",NFILE,fnm),
 835                        opt2fn_null("-ee",NFILE,fnm),
 836                        opt2fn("-m",NFILE,fnm),skip,tref,oenv);
 837
 838   if (bSplit) {
 839     printf("Splitting set of replicas in two halves\n");
 840     for(i=0; (i<remd.nreplica); i++)
 841       remd.bMask[i] = FALSE;
 842     remd.nmask = 0;
 843     for(i=0; (i<remd.nreplica); i+=2) {
 844       remd.bMask[i] = TRUE;
 845       remd.nmask++;
 846     }
 847     sum_ft(&remd);
 848     optimize_remd_parameters(fp,&remd,maxiter,tol);
 849     dump_remd_parameters(fp,&remd,"test1.xvg",NULL,NULL,NULL,NULL,skip,tref,oenv);
 850
 851     for(i=0; (i<remd.nreplica); i++)
 852       remd.bMask[i] = !remd.bMask[i];
 853     remd.nmask = remd.nreplica - remd.nmask;
 854
 855     sum_ft(&remd);
 856     optimize_remd_parameters(fp,&remd,maxiter,tol);
 857     dump_remd_parameters(fp,&remd,"test2.xvg",NULL,NULL,NULL,NULL,skip,tref,oenv);
 858
 859     for(i=0; (i<remd.nreplica); i++)
 860       remd.bMask[i] = FALSE;
 861     remd.nmask = 0;
 862     for(i=0; (i<remd.nreplica/2); i++) {
 863       remd.bMask[i] = TRUE;
 864       remd.nmask++;
 865     }
 866     sum_ft(&remd);
 867     optimize_remd_parameters(fp,&remd,maxiter,tol);
 868     dump_remd_parameters(fp,&remd,"test1.xvg",NULL,NULL,NULL,NULL,skip,tref,oenv);
 869
 870     for(i=0; (i<remd.nreplica); i++)
 871       remd.bMask[i] = FALSE;
 872     remd.nmask = 0;
 873     for(i=remd.nreplica/2; (i<remd.nreplica); i++) {
 874       remd.bMask[i] = TRUE;
 875       remd.nmask++;
 876     }
 877     sum_ft(&remd);
 878     optimize_remd_parameters(fp,&remd,maxiter,tol);
 879     dump_remd_parameters(fp,&remd,"test1.xvg",NULL,NULL,NULL,NULL,skip,tref,oenv);
 880   }
 881   ffclose(fp);
 882
 883   view_all(oenv, NFILE, fnm);
 884
 885   thanx(stderr);
 886
 887   return 0;
 888 }
 889
 890 #else
 891 int gmx_kinetics(int argc,char *argv[])
 892 {
 893   fprintf(stderr,"This program should be compiled with the GNU scientific library. Please install the library and reinstall GROMACS.\n");
 894
 895   return 0;
 896 }
 897 #endif