2 * This file is part of the GROMACS molecular simulation package.
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, 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.
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.
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.
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.
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.
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.
49 #include "gmx_fatal.h"
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" };
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. */
73 gmx_bool bDiscrete; /* Use either discrete folding (0/1) or a continuous */
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;
87 #include <gsl/gsl_multimin.h>
89 static char *itoa(int i)
97 static char *epnm(int nparams,int index)
99 static char buf[32],from[8],to[8];
102 range_check(index,0,nparams);
103 if ((nparams == 2) || (nparams == 4))
105 else if ((nparams > 4) && (nparams % 4 == 0))
108 gmx_fatal(FARGS,"Don't know how to handle %d parameters",nparams);
113 static gmx_bool bBack(t_remd_data *d)
115 return (d->nparams > 2);
118 static real is_folded(t_remd_data *d,int irep,int jframe)
120 if (d->state[irep][jframe] == 0)
126 static real is_unfolded(t_remd_data *d,int irep,int jframe)
128 if (d->state[irep][jframe] == d->nstate-1)
134 static real is_intermediate(t_remd_data *d,int irep,int jframe)
136 if ((d->state[irep][jframe] == 1) && (d->nstate > 2))
142 static void integrate_dfdt(t_remd_data *d)
145 double beta,ddf,ddi,df,db,fac,sumf,sumi,area;
149 for(i=0; (i<d->nreplica); i++) {
152 ddf = 0.5*d->dt*is_folded(d,i,0);
153 ddi = 0.5*d->dt*is_intermediate(d,i,0);
156 ddf = 0.5*d->dt*d->state[i][0];
159 d->fcalt[i][0] = ddf;
160 d->icalt[i][0] = ddi;
165 for(j=1; (j<d->nframe); j++) {
171 for(i=0; (i<d->nreplica); i++) {
173 beta = d->beta[i][j];
174 if ((d->nstate <= 2) || d->bDiscrete) {
176 df = (d->params[epAuf]*exp(-beta*d->params[epEuf])*
179 area = (d->data2 ? d->data2[i][j] : 1.0);
180 df = area*d->params[epAuf]*exp(-beta*d->params[epEuf]);
185 db = (d->params[epAfu]*exp(-beta*d->params[epEfu])*
188 gmx_fatal(FARGS,"Back reaction not implemented with continuous");
193 d->fcalt[i][j] = d->fcalt[i][j-1] + ddf;
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])*
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;
212 d->sumfct[j] = d->sumfct[j-1] + sumf;
213 d->sumict[j] = d->sumict[j-1] + sumi;
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]);
220 fprintf(debug,"&\n");
224 static void sum_ft(t_remd_data *d)
229 for(j=0; (j<d->nframe); j++) {
232 if ((j == 0) || (j==d->nframe-1))
236 for(i=0; (i<d->nreplica); i++) {
239 d->sumft[j] += fac*is_folded(d,i,j);
240 d->sumit[j] += fac*is_intermediate(d,i,j);
243 d->sumft[j] += fac*d->state[i][j];
249 static double calc_d2(t_remd_data *d)
252 double dd2,d2=0,dr2,tmp;
257 for(j=d->j0; (j<d->j1); j++) {
259 d2 += sqr(d->sumft[j]-d->sumfct[j]);
261 d2 += sqr(d->sumit[j]-d->sumict[j]);
264 d2 += sqr(d->sumft[j]-d->sumfct[j]);
268 for(i=0; (i<d->nreplica); i++) {
271 for(j=d->j0; (j<d->j1); j++) {
272 tmp = sqr(is_folded(d,i,j)-d->fcalt[i][j]);
276 tmp = sqr(is_intermediate(d,i,j)-d->icalt[i][j]);
281 d->d2_replica[i] = dr2/(d->j1-d->j0);
285 dd2 = (d2/(d->j1-d->j0))/(d->bDiscrete ? d->nmask : 1);
290 static double my_f(const gsl_vector *v,void *params)
292 t_remd_data *d = (t_remd_data *) params;
296 for(i=0; (i<d->nparams); i++) {
297 d->params[i] = gsl_vector_get(v, i);
298 if (d->params[i] < 0)
307 static void optimize_remd_parameters(FILE *fp,t_remd_data *d,int maxiter,
315 const gsl_multimin_fminimizer_type *T;
316 gsl_multimin_fminimizer *s;
319 gsl_multimin_function my_func;
322 my_func.n = d->nparams;
323 my_func.params = (void *) d;
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]);
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]);
335 T = gsl_multimin_fminimizer_nmsimplex;
336 s = gsl_multimin_fminimizer_alloc (T, my_func.n);
338 gsl_multimin_fminimizer_set (s, &my_func, x, dx);
340 gsl_vector_free (dx);
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");
349 status = gsl_multimin_fminimizer_iterate (s);
352 gmx_fatal(FARGS,"Something went wrong in the iteration in minimizer %s",
353 gsl_multimin_fminimizer_name(s));
355 d2 = gsl_multimin_fminimizer_minimum(s);
356 size = gsl_multimin_fminimizer_size(s);
357 status = gsl_multimin_test_size(size,tol);
359 if (status == GSL_SUCCESS)
360 printf ("Minimum found using %s at:\n",
361 gsl_multimin_fminimizer_name(s));
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);
368 while ((status == GSL_CONTINUE) && (iter < maxiter));
370 gsl_multimin_fminimizer_free (s);
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,
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);
390 d->nparams = 4*(1+ninter);
391 d->nstate = 2+ninter;
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);
407 /* convert_times(d->nframe,d->time); */
412 for(d->j0=0; (d->j0<d->nframe) && (d->time[d->j0] < t0); d->j0++)
417 for(d->j1=0; (d->j1<d->nframe) && (d->time[d->j1] < t1); d->j1++)
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);
430 for(j=0; (j<d->nframe); j++) {
431 d->beta[i][j] = 1.0/(BOLTZ*d->temp[i][j]);
436 else if ((ucut > cutoff) && (dd <= ucut))
439 d->state[i][j] = d->nstate-1;
442 d->state[i][j] = dd*nmult;
447 /* Assume forward rate constant is half the total time in this
448 * simulation and backward is ten times as long */
450 tau_f = d->time[d->nframe-1];
452 d->params[epEuf] = Euf;
453 d->params[epAuf] = exp(d->params[epEuf]/(BOLTZ*tref))/tau_f;
455 d->params[epEfu] = Efu;
456 d->params[epAfu] = exp(d->params[epEfu]/(BOLTZ*tref))/tau_u;
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;
465 d->params[epAfu] = 0;
466 d->params[epEfu] = 0;
470 d->params[epEuf] = Euf;
472 d->params[epAuf] = 0.1;
474 d->params[epAuf] = 20.0;
478 static real tau(real A,real E,real T)
480 return exp(E/(BOLTZ*T))/A;
483 static real folded_fraction(t_remd_data *d,real tref)
487 tauf = tau(d->params[epAuf],d->params[epEuf],tref);
488 taub = tau(d->params[epAfu],d->params[epEfu],tref);
490 return (taub/(tauf+taub));
493 static void print_tau(FILE *gp,t_remd_data *d,real tref)
495 real tauf,taub,ddd,fff,DG,DH,TDS,Tm,Tb,Te,Fb,Fe,Fm;
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]);
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];
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",
522 Fb = folded_fraction(d,Tb);
523 Fe = folded_fraction(d,Te);
524 while((Te-Tb > 0.001) && (Fm != 0.5)) {
526 Fm = folded_fraction(d,Tm);
536 if ((Fb-0.5)*(Fe-0.5) <= 0)
537 fprintf(gp,"Melting temperature Tm = %8.3f K\n",Tm);
539 fprintf(gp,"No melting temperature detected between 260 and 420K\n");
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);
553 fprintf(gp,"Equilibrium properties at T = %g\n",tref);
554 fprintf(gp,"tau_f = %8.3f\n",tauf);
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,
564 int i,j,np=d->nparams;
565 real rhs,tauf,taub,fff,DG;
567 const char *leg[] = { "Measured", "Fit", "Difference" };
568 const char *mleg[] = { "Folded fraction","DG (kJ/mole)"};
570 real fac[] = { 0.97, 0.98, 0.99, 1.0, 1.01, 1.02, 1.03 };
571 #define NFAC asize(fac)
576 print_tau(gp,d,tref);
577 norm = (d->bDiscrete ? 1.0/d->nmask : 1.0);
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);
591 if (!d->bSum && rfn) {
592 snew(rleg,d->nreplica*2);
593 for(i=0; (i<d->nreplica); i++) {
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);
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]);
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);
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);
640 snew(params,d->nparams);
641 for(i=0; (i<d->nparams); i++)
642 params[i] = d->params[i];
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];
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]);
661 for(i=0; (i<d->nparams); i++)
662 d->params[i] = params[i];
666 for(i=0; (i<d->nreplica); i++)
667 fprintf(gp,"Chi2[%3d] = %8.2e\n",i,d->d2_replica[i]);
670 #endif /*HAVE_LIBGSL*/
672 int gmx_kinetics(int argc,char *argv[])
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."
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;
705 static gmx_bool bHaveT = TRUE;
710 static real tol = 1e-3;
711 static int maxiter = 100;
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;
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" }
760 #define NPA asize(pa)
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;
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 }
780 #define NFILE asize(fnm)
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);
787 please_cite(stdout,"Spoel2006d");
789 gmx_fatal(FARGS,"cutoff should be >= 0 (rather than %f)",cutoff);
791 tfile = opt2fn("-f",NFILE,fnm);
792 dfile = opt2fn("-d",NFILE,fnm);
793 dfile2 = opt2fn_null("-d2",NFILE,fnm);
795 fp = ffopen(opt2fn("-g",NFILE,fnm),"w");
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);
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);
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",
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",
829 remd.nreplica = nset_d;
832 preprocess_remd(fp,&remd,cutoff,tref,ucut,bBack,Euf,Efu,Ei,t0,t1,
833 bSum,bDiscrete,nmult);
835 optimize_remd_parameters(fp,&remd,maxiter,tol);
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);
844 printf("Splitting set of replicas in two halves\n");
845 for(i=0; (i<remd.nreplica); i++)
846 remd.bMask[i] = FALSE;
848 for(i=0; (i<remd.nreplica); i+=2) {
849 remd.bMask[i] = TRUE;
853 optimize_remd_parameters(fp,&remd,maxiter,tol);
854 dump_remd_parameters(fp,&remd,"test1.xvg",NULL,NULL,NULL,NULL,skip,tref,oenv);
856 for(i=0; (i<remd.nreplica); i++)
857 remd.bMask[i] = !remd.bMask[i];
858 remd.nmask = remd.nreplica - remd.nmask;
861 optimize_remd_parameters(fp,&remd,maxiter,tol);
862 dump_remd_parameters(fp,&remd,"test2.xvg",NULL,NULL,NULL,NULL,skip,tref,oenv);
864 for(i=0; (i<remd.nreplica); i++)
865 remd.bMask[i] = FALSE;
867 for(i=0; (i<remd.nreplica/2); i++) {
868 remd.bMask[i] = TRUE;
872 optimize_remd_parameters(fp,&remd,maxiter,tol);
873 dump_remd_parameters(fp,&remd,"test1.xvg",NULL,NULL,NULL,NULL,skip,tref,oenv);
875 for(i=0; (i<remd.nreplica); i++)
876 remd.bMask[i] = FALSE;
878 for(i=remd.nreplica/2; (i<remd.nreplica); i++) {
879 remd.bMask[i] = TRUE;
883 optimize_remd_parameters(fp,&remd,maxiter,tol);
884 dump_remd_parameters(fp,&remd,"test1.xvg",NULL,NULL,NULL,NULL,skip,tref,oenv);
888 view_all(oenv, NFILE, fnm);
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*/