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