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