From: Roland Schulz <roland@utk.edu>
Date: Mon, 1 Oct 2012 13:20:00 +0000 (-0400)
Subject: Merge release-4-5-patches into release-4-6
X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=597c9d5a798e0141d5eeaf602a2069b21e335d1a;p=alexxy%2Fgromacs.git

Merge release-4-5-patches into release-4-6

Conflicts:
	src/gmxlib/confio.c
	src/gmxlib/copyrite.c
	src/gmxlib/statutil.c
	src/kernel/readir.c
	src/kernel/repl_ex.c
	src/kernel/runner.c

Change-Id: I3828c50f8166a6095ee64672cfd4380d3e4c86f6
---

597c9d5a798e0141d5eeaf602a2069b21e335d1a
diff --cc src/gmxlib/confio.c
index a1112cf40c,41a509e307..886983a6d5
--- a/src/gmxlib/confio.c
+++ b/src/gmxlib/confio.c
@@@ -880,9 -878,10 +878,10 @@@ static void read_whole_conf(const char 
  gmx_bool gro_next_x_or_v(FILE *status,t_trxframe *fr)
  {
    t_atoms atoms;
+   t_symtab symtab;
    char    title[STRLEN],*p;
    double  tt;
 -  int     ndec,i;
 +  int     ndec=0,i;
  
    if (gmx_eof(status))
      return FALSE;
diff --cc src/gmxlib/copyrite.c
index 5e1d0a6e29,d01a9c0a1b..41bd02830c
--- a/src/gmxlib/copyrite.c
+++ b/src/gmxlib/copyrite.c
@@@ -537,32 -527,13 +537,37 @@@ void please_cite(FILE *fp,const char *k
        "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)
    
diff --cc src/gmxlib/statutil.c
index f1d79985e2,7554102982..fd7ff31cae
--- a/src/gmxlib/statutil.c
+++ b/src/gmxlib/statutil.c
@@@ -252,8 -252,8 +252,8 @@@ int check_times(real t
  
  static void set_default_time_unit(const char *time_list[], gmx_bool bCanTime)
  {
 -    int i,j;
 +    int i=0,j;
-     const char *select;
+     const char *select = NULL;
  
      if (bCanTime)
      {
diff --cc src/kernel/readir.c
index ff54f6b3ce,8163d9275b..24451eaaa8
--- a/src/kernel/readir.c
+++ b/src/kernel/readir.c
@@@ -710,7 -431,10 +710,7 @@@ void check_ir(const char *mdparin,t_inp
                (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 &&
                 ir->compress[ZZ][XX] <= 0 && ir->compress[ZZ][YY] <= 0));
          
-         if (epcPARRINELLORAHMAN == ir->epct && opts->bGenVel)
 -        sprintf(err_buf,"pressure coupling with PPPM not implemented, use PME");
 -        CHECK(ir->coulombtype == eelPPPM);
 -
+         if (epcPARRINELLORAHMAN == ir->epc && opts->bGenVel)
          {
              sprintf(warn_buf,
                      "You are generating velocities so I am assuming you "
diff --cc src/kernel/repl_ex.c
index d20b189df3,1ee8998d21..8b6797051b
--- a/src/kernel/repl_ex.c
+++ b/src/kernel/repl_ex.c
@@@ -171,13 -176,11 +171,13 @@@ gmx_repl_ex_t init_replica_exchange(FIL
      {
          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;
@@@ -727,148 -584,22 +727,148 @@@ static real calc_delta(FILE *fplog, gmx
      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);
      }
 +    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,
-                                  int step,real time,int *pind)
++                                 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 %d time %g\n",step,time);
++    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);
  
@@@ -1030,152 -716,24 +1030,152 @@@ static void write_debug_x(t_state *stat
      }
  }
  
 +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
diff --cc src/kernel/repl_ex.h
index bc5006d13d,0d81016461..7e4bf23d68
--- a/src/kernel/repl_ex.h
+++ b/src/kernel/repl_ex.h
@@@ -51,9 -51,9 +51,9 @@@ extern gmx_repl_ex_t init_replica_excha
  extern gmx_bool replica_exchange(FILE *fplog,
  			     const t_commrec *cr,
  			     gmx_repl_ex_t re,
 -			     t_state *state,real *ener,
 +			     t_state *state,gmx_enerdata_t *enerd,
  			     t_state *state_local,
- 			     int step,real time);
+ 			     gmx_large_int_t step,real time);
  /* Attempts replica exchange, should be called on all nodes.
   * Returns TRUE if this state has been exchanged.
   * When running each replica in parallel,
diff --cc src/mdlib/gmx_wallcycle.c
index 8188c177f7,3fa71b84fe..b62f921c2a
--- a/src/mdlib/gmx_wallcycle.c
+++ b/src/mdlib/gmx_wallcycle.c
@@@ -436,10 -430,10 +436,10 @@@ void wallcycle_print(FILE *fplog, int n
          fprintf(fplog,"%s\n",myline);
          for(i=ewcPPDURINGPME+1; i<ewcNR; i++)
          {
 -            if (subdivision(i))
 +            if (pme_subdivision(i))
              {
                  print_cycles(fplog,c2t,wcn[i],
-                              (i>=ewcPMEMESH || i<=ewcPME_SOLVE) ? npme : npp,
+                              (i>=ewcPMEMESH && i<=ewcPME_SOLVE) ? npme : npp,
                               wc->wcc[i].n,cycles[i],tot);
              }
          }