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