"H. Wang, F. Dommert, C.Holm",
"Optimizing working parameters of the smooth particle mesh Ewald algorithm in terms of accuracy and efficiency",
"J. Chem. Phys. B",
- 133, 2010, "034117"
- },
+ 133, 2010, "034117" },
+ { "Sugita1999a",
+ "Y. Sugita, Y. Okamoto",
+ "Replica-exchange molecular dynamics method for protein folding",
+ "Chem. Phys. Lett.",
+ 314, 1999, "141-151" },
+ { "Kutzner2011",
+ "C. Kutzner and J. Czub and H. Grubmuller",
+ "Keep it Flexible: Driving Macromolecular Rotary Motions in Atomistic Simulations with GROMACS",
+ "J. Chem. Theory Comput.",
+ 7, 2011, "1381-1393" },
+ { "Hoefling2011",
+ "M. Hoefling, N. Lima, D. Haenni, C.A.M. Seidel, B. Schuler, H. Grubmuller",
+ "Structural Heterogeneity and Quantitative FRET Efficiency Distributions of Polyprolines through a Hybrid Atomistic Simulation and Monte Carlo Approach",
+ "PLoS ONE",
+ 6, 2011, "e19791" },
+ { "Hockney1988",
+ "R. W. Hockney and J. W. Eastwood",
+ "Computer simulation using particles",
+ "IOP, Bristol",
+ 1, 1988, "1" },
+ { "Ballenegger2012",
+ "V. Ballenegger, J.J. Cerda, and C. Holm",
+ "How to Convert SPME to P3M: Influence Functions and Error Estimates",
+ "J. Chem. Theory Comput.",
+ 8, 2012, "936-947" },
+ { "Garmay2012",
+ "Garmay Yu, Shvetsov A, Karelov D, Lebedev D, Radulescu A, Petukhov M, Isaev-Ivanov V",
+ "Correlated motion of protein subdomains and large-scale conformational flexibility of RecA protein filament",
+ "Journal of Physics: Conference Series",
+ 340, 2012, "012094" }
};
#define NSTR (int)asize(citedb)
{
gmx_fatal(FARGS,"The properties of the %d systems are all the same, there is nothing to exchange",re->nrepl);
}
+ if (bLambda && bTemp) {
+ re->type = ereTL;
+ }
- switch (re->type)
+ if (bTemp)
{
- please_cite(fplog,"Hukushima96a");
- case ereTEMP:
+ please_cite(fplog,"Sugita1999a");
if (ir->epc != epcNO)
{
re->bNPT = TRUE;
switch (re->type)
{
case ereTEMP:
- snew(Epot,re->nrepl);
- snew(Vol,re->nrepl);
- Epot[re->repl] = ener[F_EPOT];
- Vol[re->repl] = vol;
- gmx_sum_sim(re->nrepl,Epot,ms);
- gmx_sum_sim(re->nrepl,Vol,ms);
+ /*
+ * Okabe et. al. Chem. Phys. Lett. 335 (2001) 435-439
+ */
+ ediff = Epot[b] - Epot[a];
+ delta = -(beta[bp] - beta[ap])*ediff;
break;
case ereLAMBDA:
- snew(dvdl,re->nrepl);
- dvdl[re->repl] = ener[F_DVDL];
- gmx_sum_sim(re->nrepl,dvdl,ms);
+ /* two cases: when we are permuted, and not. */
+ /* non-permuted:
+ ediff = E_new - E_old
+ = [H_b(x_a) + H_a(x_b)] - [H_b(x_b) + H_a(x_a)]
+ = [H_b(x_a) - H_a(x_a)] + [H_a(x_b) - H_b(x_b)]
+ = df[b][a] + df[a][b] */
+ /* permuted:
+ ediff = E_new - E_old
+ = [H_bp(x_a) + H_ap(x_b)] - [H_bp(x_b) + H_ap(x_a)]
+ = [H_bp(x_a) - H_ap(x_a)] + [H_ap(x_b) - H_bp(x_b)]
+ = [H_bp(x_a) - H_a(x_a) + H_a(x_a) - H_ap(x_a)] + [H_ap(x_b) - H_b(x_b) + H_b(x_b) - H_bp(x_b)]
+ = [H_bp(x_a) - H_a(x_a)] - [H_ap(x_a) - H_a(x_a)] + [H_ap(x_b) - H_b(x_b)] - H_bp(x_b) - H_b(x_b)]
+ = (df[bp][a] - df[ap][a]) + (df[ap][b] - df[bp][b]) */
+ ediff = (df[bp][a] - df[ap][a]) + (df[ap][b] - df[bp][b]);
+ delta = ediff*beta[a]; /* assume all same temperature in this case */
break;
+ case ereTL:
+ /* not permuted: */
+ /* delta = reduced E_new - reduced E_old
+ = [beta_b H_b(x_a) + beta_a H_a(x_b)] - [beta_b H_b(x_b) + beta_a H_a(x_a)]
+ = [beta_b H_b(x_a) - beta_a H_a(x_a)] + [beta_a H_a(x_b) - beta_b H_b(x_b)]
+ = [beta_b dH_b(x_a) + beta_b H_a(x_a) - beta_a H_a(x_a)] +
+ [beta_a dH_a(x_b) + beta_a H_b(x_b) - beta_b H_b(x_b)]
+ = [beta_b dH_b(x_a) + [beta_a dH_a(x_b) +
+ beta_b (H_a(x_a) - H_b(x_b)]) - beta_a (H_a(x_a) - H_b(x_b))
+ = beta_b dH_b(x_a) + beta_a dH_a(x_b) - (beta_b - beta_a)(H_b(x_b) - H_a(x_a) */
+ /* delta = beta[b]*df[b][a] + beta[a]*df[a][b] - (beta[b] - beta[a])*(Epot[b] - Epot[a]; */
+ /* permuted (big breath!) */
+ /* delta = reduced E_new - reduced E_old
+ = [beta_bp H_bp(x_a) + beta_ap H_ap(x_b)] - [beta_bp H_bp(x_b) + beta_ap H_ap(x_a)]
+ = [beta_bp H_bp(x_a) - beta_ap H_ap(x_a)] + [beta_ap H_ap(x_b) - beta_bp H_bp(x_b)]
+ = [beta_bp H_bp(x_a) - beta_ap H_ap(x_a)] + [beta_ap H_ap(x_b) - beta_bp H_bp(x_b)]
+ - beta_pb H_a(x_a) + beta_ap H_a(x_a) + beta_pb H_a(x_a) - beta_ap H_a(x_a)
+ - beta_ap H_b(x_b) + beta_bp H_b(x_b) + beta_ap H_b(x_b) - beta_bp H_b(x_b)
+ = [(beta_bp H_bp(x_a) - beta_bp H_a(x_a)) - (beta_ap H_ap(x_a) - beta_ap H_a(x_a))] +
+ [(beta_ap H_ap(x_b) - beta_ap H_b(x_b)) - (beta_bp H_bp(x_b) - beta_bp H_b(x_b))]
+ + beta_pb H_a(x_a) - beta_ap H_a(x_a) + beta_ap H_b(x_b) - beta_bp H_b(x_b)
+ = [beta_bp (H_bp(x_a) - H_a(x_a)) - beta_ap (H_ap(x_a) - H_a(x_a))] +
+ [beta_ap (H_ap(x_b) - H_b(x_b)) - beta_bp (H_bp(x_b) - H_b(x_b))]
+ + beta_pb (H_a(x_a) - H_b(x_b)) - beta_ap (H_a(x_a) - H_b(x_b))
+ = ([beta_bp df[bp][a] - beta_ap df[ap][a]) + beta_ap df[ap][b] - beta_bp df[bp][b])
+ + (beta_pb-beta_ap)(H_a(x_a) - H_b(x_b)) */
+ delta = beta[bp]*(df[bp][a] - df[bp][b]) + beta[ap]*(df[ap][b] - df[ap][a]) - (beta[bp]-beta[ap])*(Epot[b]-Epot[a]);
+ break;
+ default:
+ gmx_incons("Unknown replica exchange quantity");
+ }
+ if (bPrint)
+ {
+ fprintf(fplog,"Repl %d <-> %d dE_term = %10.3e (kT)\n",a,b,delta);
}
- int step,real time,int *pind)
+ if (re->bNPT)
+ {
+ /* revist the calculation for 5.0. Might be some improvements. */
+ dpV = (beta[ap]*re->pres[ap]-beta[bp]*re->pres[bp])*(Vol[b]-Vol[a])/PRESFAC;
+ if (bPrint)
+ {
+ fprintf(fplog," dpV = %10.3e d = %10.3e\nb",dpV,delta + dpV);
+ }
+ delta += dpV;
+ }
+ return delta;
+}
+
+static void get_replica_exchange(FILE *fplog,const gmx_multisim_t *ms,
+ struct gmx_repl_ex *re,gmx_enerdata_t *enerd,real vol,
- fprintf(fplog,"Replica exchange at step %d time %g\n",step,time);
++ gmx_large_int_t step,real time,int *pind)
+{
+ int m,i,a,b,ap,bp,i0,i1,tmp;
+ real *Epot=NULL,*Vol=NULL,**flambda=NULL,*beta=NULL,*prob;
+ real ediff=0,delta=0,dpV=0;
+ gmx_bool *bEx,bPrint,bMultiEx;
+ gmx_bool bEpot=FALSE;
+ gmx_bool bFLambda=FALSE;
+ gmx_bool bVol=FALSE;
+
+ bMultiEx = (re->nex > 1); /* multiple exchanges at each state */
++ fprintf(fplog,"Replica exchange at step " gmx_large_int_pfmt " time %g\n",step,time);
+ snew(beta,re->nrepl);
+ if (re->bNPT)
+ {
+ bVol = TRUE;
+ snew(Vol,re->nrepl);
+ Vol[re->repl] = vol;
+ }
+
+ if ((re->type == ereTEMP || re->type == ereTL))
+ {
+ bEpot = TRUE;
+ snew(Epot,re->nrepl);
+ Epot[re->repl] = enerd->term[F_EPOT];
+ /* temperatures of different states*/
+ for (i=0;i<re->nrepl;i++)
+ {
+ beta[i] = 1.0/(re->q[ereTEMP][i]*BOLTZ);
+ }
+ }
+ else
+ {
+ for (i=0;i<re->nrepl;i++)
+ {
+ beta[i] = 1.0/(re->temp*BOLTZ); /* we have a single temperature */
+ }
+ }
+ if (re->type == ereLAMBDA || re->type == ereTL)
+ {
+ bFLambda = TRUE;
+ /* lambda differences. */
+ /* flambda[i][j] is the energy of the jth simulation in the ith Hamiltonian
+ minus the energy of the jth simulation in the jth Hamiltonian */
+ snew(flambda,re->nrepl);
+ for (i=0;i<re->nrepl;i++)
+ {
+ snew(flambda[i],re->nrepl);
+ flambda[i][re->repl] = (enerd->enerpart_lambda[(int)re->q[ereLAMBDA][i]+1]-enerd->enerpart_lambda[0]);
+ }
+ }
+
+ /* now actually do the communication */
+ if (bVol)
+ {
+ gmx_sum_sim(re->nrepl,Vol,ms);
+ }
+ if (bEpot)
+ {
+ gmx_sum_sim(re->nrepl,Epot,ms);
+ }
+ if (bFLambda)
+ {
+ for (i=0;i<re->nrepl;i++)
+ {
+ gmx_sum_sim(re->nrepl,flambda[i],ms);
+ }
+ }
snew(bEx,re->nrepl);
snew(prob,re->nrepl);
}
}
+static void cyclic_decomposition(FILE *fplog, int *pind, int **cyclic, int n, int *nswap)
+{
+
+ int i,j,c,p;
+ int *incycle;
+ int maxlen = 1;
+ snew(incycle,n);
+
+ /* compute cyclic decompositions */
+ for (i=0;i<n;i++)
+ {
+ snew(cyclic[i],n);
+ for (j=0;j<n;j++)
+ {
+ cyclic[i][j] = -1;
+ }
+ }
+
+ for (i=0;i<n;i++) /* one cycle for each replica */
+ {
+ if (incycle[i])
+ {
+ cyclic[i][0] = -1;
+ continue;
+ }
+ cyclic[i][0] = i;
+ incycle[i] = TRUE;
+ c = 1;
+ p = i;
+ for (j=0;j<n;j++) /* potentially all cycles are part, but we will break first */
+ {
+ p = pind[p]; /* start permuting */
+ if (p==i)
+ {
+ cyclic[i][c] = -1;
+ if (c > maxlen)
+ {
+ maxlen = c;
+ }
+ break; /* we've reached the original element, the cycle is complete, and we marked the end. */
+ }
+ else
+ {
+ cyclic[i][c] = p; /* each permutation gives a new member of the cycle */
+ incycle[p] = TRUE;
+ c++;
+ }
+ }
+ }
+ *nswap = maxlen - 1;
+
+ if (debug)
+ {
+ for (i=0;i<n;i++)
+ {
+ fprintf(fplog,"Cycle %d:",i);
+ for (j=0;j<n;j++)
+ {
+ if (cyclic[i][j] < 0)
+ {
+ break;
+ }
+ fprintf(fplog,"%2d",cyclic[i][j]);
+ }
+ fprintf(fplog,"\n");
+ }
+ fflush(fplog);
+ }
+}
+
+static void compute_exchange_order(FILE *fplog, int **cyclic,int **order, int n, int maxswap)
+{
+ int i,j;
+
+ for (i=0;i<n;i++)
+ {
+ snew(order[i],maxswap);
+ for (j=0;j<maxswap;j++)
+ {
+ order[i][j] = -1;
+ }
+ }
+ for (j=0;j<maxswap;j++)
+ {
+ for (i=0;i<n;i++)
+ {
+ if (cyclic[i][j+1] >= 0)
+ {
+ order[cyclic[i][j+1]][j] = cyclic[i][j];
+ order[cyclic[i][j]][j] = cyclic[i][j+1];
+ }
+ }
+ for (i=0;i<n;i++)
+ {
+ if (order[i][j] < 0)
+ {
+ order[i][j] = i; /* if it's not exchanging, it should stay this round*/
+ }
+ }
+ }
+ if (debug)
+ {
+ fprintf(fplog,"Replica Exchange Order\n");
+ for (i=0;i<n;i++)
+ {
+ fprintf(fplog,"Replica %d:",i);
+ for (j=0;j<maxswap;j++)
+ {
+ if (order[i][j] < 0) break;
+ fprintf(fplog,"%2d",order[i][j]);
+ }
+ fprintf(fplog,"\n");
+ }
+ fflush(fplog);
+ }
+}
+
gmx_bool replica_exchange(FILE *fplog,const t_commrec *cr,struct gmx_repl_ex *re,
- t_state *state,real *ener,
- t_state *state_local,
- gmx_large_int_t step,real time)
+ t_state *state,gmx_enerdata_t *enerd,
- t_state *state_local,int step,real time)
++ t_state *state_local,gmx_large_int_t step,real time)
{
gmx_multisim_t *ms;
- int exchange=-1,shift;
+ int exchange=-1,shift;
+ int i,j,maxswap=0;
+ int *exchanges=NULL;
+ int **cyclic=NULL;
+ int **order=NULL;
gmx_bool bExchanged=FALSE;
-
+
ms = cr->ms;
-
if (MASTER(cr))
{
- exchange = get_replica_exchange(fplog,ms,re,ener,det(state->box),
- step,time);
- bExchanged = (exchange >= 0);
+ snew(exchanges,re->nrepl);
+ get_replica_exchange(fplog,ms,re,enerd,det(state->box),step,time,exchanges);
+ bExchanged = (exchanges[re->repl] != re->nrepl); /* only mark as exchanged if it has a shuffled index */
+ snew(cyclic,re->nrepl);
+ snew(order,re->nrepl);
+
+ /* identify the cyclic decomposition of the permutation (very easy if neighbor replica exchange) */
+ cyclic_decomposition(fplog,exchanges,cyclic,re->nrepl,&maxswap);
+
+ /* now translate the decompsition into a replica exchange order at each step */
+ compute_exchange_order(fplog,cyclic,order,re->nrepl,maxswap);
+
+ sfree(cyclic); /* don't need this anymore */
}
-
if (PAR(cr))
{
#ifdef GMX_MPI
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff