?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Parse stderr of tuning runs to catch fatal errors which do not appear in md.log
[gromacs.git] / src / tools / gmx_kinetics.c
blobb498f803e872a2d0314482fa84a6827242ef868f
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 bool *bMask; /* Determine whether this replica is part of the d2 comp. */
72 bool bSum;
73 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 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,bool bBack,real Euf,real Efu,
372 real Ei,real t0,real t1,bool bSum,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 char *leg[] = { "Measured", "Fit", "Difference" };
565 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,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 "g_kinetics reads two xvg 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 xvg 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 xvg 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 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 bool bBack = TRUE;
711 static bool bSplit = TRUE;
712 static bool bSum = TRUE;
713 static 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 xvg 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 gmx_fatal(FARGS,"This program should be compiled with the GNU scientific library. Please install the library and reinstall GROMACS.");
894
895 return 0;
896 }
897 #endif
The ultimate molecular dynamics simulation package