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

Change-Id: Id4fc52b784ac55c2e57669af391b8d2a70a56695

diff --git a/cmake/gmxCFlags.cmake b/cmake/gmxCFlags.cmake
index 3b9489449f98a2da771cebfed1f67ff12f3c55d2..f6dbdb126520f552813533a4ed0ea522c0e2997a 100644 (file)
--- a/cmake/gmxCFlags.cmake
+++ b/cmake/gmxCFlags.cmake
@@ -134,9 +134,9 @@ MACRO(gmx_c_flags)
         GMX_TEST_CFLAG(CFLAGS_WARN "-Wall -Wno-unused -Wunused-value" GMXC_CFLAGS)
     endif()
 
-    if (CMAKE_C_COMPILER_ID MATCHES "Clang")
+    if (CMAKE_CXX_COMPILER_ID MATCHES "Clang")
         if(NOT GMX_OPENMP)
-            GMX_TEST_CFLAG(CXXFLAGS_PRAGMA "-Wno-unknown-pragmas" GMXC_CXXFLAGS)
+            GMX_TEST_CXXFLAG(CXXFLAGS_PRAGMA "-Wno-unknown-pragmas" GMXC_CXXFLAGS)
         endif()
         GMX_TEST_CXXFLAG(CXXFLAGS_WARN "-Wall -Wno-unused -Wunused-value" GMXC_CXXFLAGS)
     endif()

diff --git a/include/vsite.h b/include/vsite.h
index bb34b26712f53fe1e4e3577e855c36b06a876082..024b6af4fbc7de9206ca69d31df84a989143943b 100644 (file)
--- a/include/vsite.h
+++ b/include/vsite.h
@@ -90,12 +90,17 @@ void construct_vsites_mtop(FILE *log,gmx_vsite_t *vsite,
  */
 
 void spread_vsite_f(FILE *log,gmx_vsite_t *vsite,
-                          rvec x[],rvec f[],rvec *fshift,
-                          t_nrnb *nrnb,t_idef *idef,
-                          int ePBC,gmx_bool bMolPBC,t_graph *g,matrix box,
-                          t_commrec *cr);
+                   rvec x[],rvec f[],rvec *fshift,
+                   gmx_bool VirCorr,matrix vir,
+                   t_nrnb *nrnb,t_idef *idef,
+                   int ePBC,gmx_bool bMolPBC,t_graph *g,matrix box,
+                   t_commrec *cr);
 /* Spread the force operating on the vsite atoms on the surrounding atoms.
  * If fshift!=NULL also update the shift forces.
+ * If VirCorr=TRUE add the virial correction for non-linear vsite constructs
+ * to vir. This correction is required when the virial is not calculated
+ * afterwards from the particle position and forces, but in a different way,
+ * as for instance for the PME mesh contribution.
  */
 
 gmx_vsite_t *init_vsite(gmx_mtop_t *mtop,t_commrec *cr);

diff --git a/src/kernel/readir.c b/src/kernel/readir.c
index ad269e35a2d1d31778fc0d4df37e0b1189b7a142..fe55fc76fc7b5445e5ebae2b30c09932885d9459 100644 (file)
--- a/src/kernel/readir.c
+++ b/src/kernel/readir.c
@@ -434,7 +434,19 @@ void check_ir(const char *mdparin,t_inputrec *ir, t_gromppopts *opts,
         
         sprintf(err_buf,"pressure coupling with PPPM not implemented, use PME");
         CHECK(ir->coulombtype == eelPPPM);
-        
+
+        if (epcPARRINELLORAHMAN == ir->epct && opts->bGenVel)
+        {
+            sprintf(warn_buf,
+                    "You are generating velocities so I am assuming you "
+                    "are equilibrating a system. You are using "
+                    "Parrinello-Rahman pressure coupling, but this can be "
+                    "unstable for equilibration. If your system crashes, try "
+                    "equilibrating first with Berendsen pressure coupling. If "
+                    "you are not equilibrating the system, you can probably "
+                    "ignore this warning.");
+            warning(wi,warn_buf);
+        }
     }
     else if (ir->coulombtype == eelPPPM)
     {

diff --git a/src/mdlib/constr.c b/src/mdlib/constr.c
index 815076e60f91221c8f6c10fa9fdefe2a84e64eac..94abc70af7d517fe6330ee33dcf85caa4c5d86a5 100644 (file)
--- a/src/mdlib/constr.c
+++ b/src/mdlib/constr.c
@@ -773,8 +773,6 @@ void set_constraints(struct gmx_constr *constr,
                 constr->lagr_nalloc = over_alloc_dd(ncons);
                 srenew(constr->lagr,constr->lagr_nalloc);
             }
-
-            constr->shaked = shake_init();
         }
     }
 
@@ -1035,6 +1033,8 @@ gmx_constr_t init_constraints(FILE *fplog,
             {
                 please_cite(fplog,"Barth95a");
             }
+
+            constr->shaked = shake_init();
         }
     }
   

diff --git a/src/mdlib/coupling.c b/src/mdlib/coupling.c
index feb2bdeb47cdddfb6619b03b5af0af12079ce158..0e822326b89f2aaf18486b84ce167e7bbd4203a1 100644 (file)
--- a/src/mdlib/coupling.c
+++ b/src/mdlib/coupling.c
@@ -417,7 +417,12 @@ void parrinellorahman_pcoupl(FILE *fplog,gmx_large_int_t step,
     
     if (maxchange > 0.01 && fplog) {
       char buf[22];
-      fprintf(fplog,"\nStep %s  Warning: Pressure scaling more than 1%%.\n",
+      fprintf(fplog,
+              "\nStep %s  Warning: Pressure scaling more than 1%%. "
+              "This may mean your system\n is not yet equilibrated. "
+              "Use of Parrinello-Rahman pressure coupling during\n"
+              "equilibration can lead to simulation instability, "
+              "and is discouraged.\n",
              gmx_step_str(step,buf));
     }
   }

diff --git a/src/mdlib/sim_util.c b/src/mdlib/sim_util.c
index 57523b5ba62519402e9e5dc10beff5da6c811f34..f5cbedd5f9380961b0fae4b40dc588d30b5ca4cc 100644 (file)
--- a/src/mdlib/sim_util.c
+++ b/src/mdlib/sim_util.c
@@ -851,14 +851,14 @@ void do_force(FILE *fplog,t_commrec *cr,
         if (vsite && !(fr->bF_NoVirSum && !(flags & GMX_FORCE_VIRIAL)))
         {
             wallcycle_start(wcycle,ewcVSITESPREAD);
-            spread_vsite_f(fplog,vsite,x,f,fr->fshift,nrnb,
+            spread_vsite_f(fplog,vsite,x,f,fr->fshift,FALSE,NULL,nrnb,
                            &top->idef,fr->ePBC,fr->bMolPBC,graph,box,cr);
             wallcycle_stop(wcycle,ewcVSITESPREAD);
 
             if (bSepLRF)
             {
                 wallcycle_start(wcycle,ewcVSITESPREAD);
-                spread_vsite_f(fplog,vsite,x,fr->f_twin,NULL,
+                spread_vsite_f(fplog,vsite,x,fr->f_twin,NULL,FALSE,NULL,
                                nrnb,
                                &top->idef,fr->ePBC,fr->bMolPBC,graph,box,cr);
                 wallcycle_stop(wcycle,ewcVSITESPREAD);
@@ -935,7 +935,9 @@ void do_force(FILE *fplog,t_commrec *cr,
              * if the constructing atoms aren't local.
              */
             wallcycle_start(wcycle,ewcVSITESPREAD);
-            spread_vsite_f(fplog,vsite,x,fr->f_novirsum,NULL,nrnb,
+            spread_vsite_f(fplog,vsite,x,fr->f_novirsum,NULL,
+                           (flags & GMX_FORCE_VIRIAL),fr->vir_el_recip,
+                           nrnb,
                            &top->idef,fr->ePBC,fr->bMolPBC,graph,box,cr);
             wallcycle_stop(wcycle,ewcVSITESPREAD);
         }

diff --git a/src/mdlib/vsite.c b/src/mdlib/vsite.c
index b2ce30d94d868429e9e307222feb7aa9bc1d651a..34a340b74f0d33d8731327670b2625bb4d7bb3c9 100644 (file)
--- a/src/mdlib/vsite.c
+++ b/src/mdlib/vsite.c
@@ -1,4 +1,5 @@
-/*
+/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
+ *
  * 
  *                This source code is part of
  * 
@@ -694,8 +695,9 @@ static void spread_vsite3(t_iatom ia[],real a,real b,
 }
 
 static void spread_vsite3FD(t_iatom ia[],real a,real b,
-                           rvec x[],rvec f[],rvec fshift[],
-                           t_pbc *pbc,t_graph *g)
+                            rvec x[],rvec f[],rvec fshift[],
+                            gmx_bool VirCorr,matrix dxdf,
+                            t_pbc *pbc,t_graph *g)
 {
   real fx,fy,fz,c,invl,fproj,a1;
   rvec xvi,xij,xjk,xix,fv,temp;
@@ -771,12 +773,38 @@ static void spread_vsite3FD(t_iatom ia[],real a,real b,
     fshift[    skj][ZZ] += a*temp[ZZ];
   }
 
+    if (VirCorr)
+    {
+        /* When VirCorr=TRUE, the virial for the current forces is not
+         * calculated from the redistributed forces. This means that
+         * the effect of non-linear virtual site constructions on the virial
+         * needs to be added separately. This contribution can be calculated
+         * in many ways, but the simplest and cheapest way is to use
+         * the first constructing atom ai as a reference position in space:
+         * subtract (xv-xi)*fv and add (xj-xi)*fj + (xk-xi)*fk.
+         */
+        rvec xiv;
+        int  i,j;
+
+        pbc_rvec_sub(pbc,x[av],x[ai],xiv);
+
+        for(i=0; i<DIM; i++)
+        {
+            for(j=0; j<DIM; j++)
+            {
+                /* As xix is a linear combination of j and k, use that here */
+                dxdf[i][j] += -xiv[i]*fv[j] + xix[i]*temp[j];
+            }
+        }
+    }
+
   /* TOTAL: 61 flops */
 }
 
 static void spread_vsite3FAD(t_iatom ia[],real a,real b,
-                            rvec x[],rvec f[],rvec fshift[],
-                            t_pbc *pbc,t_graph *g)
+                             rvec x[],rvec f[],rvec fshift[],
+                             gmx_bool VirCorr,matrix dxdf,
+                             t_pbc *pbc,t_graph *g)
 {
   rvec    xvi,xij,xjk,xperp,Fpij,Fppp,fv,f1,f2,f3;
   real    a1,b1,c1,c2,invdij,invdij2,invdp,fproj;
@@ -860,13 +888,34 @@ static void spread_vsite3FAD(t_iatom ia[],real a,real b,
     fshift[    skj][YY] +=                          f2[YY];
     fshift[    skj][ZZ] +=                          f2[ZZ];
   }
+
+    if (VirCorr)
+    {
+        rvec xiv;
+        int  i,j;
+
+        pbc_rvec_sub(pbc,x[av],x[ai],xiv);
+
+        for(i=0; i<DIM; i++)
+        {
+            for(j=0; j<DIM; j++)
+            {
+                /* Note that xik=xij+xjk, so we have to add xij*f2 */
+                dxdf[i][j] +=
+                    - xiv[i]*fv[j]
+                    + xij[i]*(f1[j] + (1 - c2)*f2[j] - f3[j])
+                    + xjk[i]*f2[j];
+            }
+        }
+    }
   
   /* TOTAL: 113 flops */
 }
 
 static void spread_vsite3OUT(t_iatom ia[],real a,real b,real c,
-                            rvec x[],rvec f[],rvec fshift[],
-                            t_pbc *pbc,t_graph *g)
+                             rvec x[],rvec f[],rvec fshift[],
+                             gmx_bool VirCorr,matrix dxdf,
+                             t_pbc *pbc,t_graph *g)
 {
   rvec    xvi,xij,xik,fv,fj,fk;
   real    cfx,cfy,cfz;
@@ -927,12 +976,29 @@ static void spread_vsite3OUT(t_iatom ia[],real a,real b,real c,
     rvec_inc(fshift[sji],fj);
     rvec_inc(fshift[ski],fk);
   }
+
+    if (VirCorr)
+    {
+        rvec xiv;
+        int  i,j;
+
+        pbc_rvec_sub(pbc,x[av],x[ai],xiv);
+
+        for(i=0; i<DIM; i++)
+        {
+            for(j=0; j<DIM; j++)
+            {
+                dxdf[i][j] += -xiv[i]*fv[j] + xij[i]*fj[j] + xik[i]*fk[j];
+            }
+        }
+    }
   
   /* TOTAL: 54 flops */
 }
 
 static void spread_vsite4FD(t_iatom ia[],real a,real b,real c,
                             rvec x[],rvec f[],rvec fshift[],
+                            gmx_bool VirCorr,matrix dxdf,
                             t_pbc *pbc,t_graph *g)
 {
   real    d,invl,fproj,a1;
@@ -1007,12 +1073,29 @@ static void spread_vsite4FD(t_iatom ia[],real a,real b,real c,
     }
   }
 
+    if (VirCorr)
+    {
+        rvec xiv;
+        int  i,j;
+
+        pbc_rvec_sub(pbc,x[av],x[ai],xiv);
+
+        for(i=0; i<DIM; i++)
+        {
+            for(j=0; j<DIM; j++)
+            {
+                dxdf[i][j] += -xiv[i]*fv[j] + xix[i]*temp[j];
+            }
+        }
+    }
+
   /* TOTAL: 77 flops */
 }
 
 
 static void spread_vsite4FDN(t_iatom ia[],real a,real b,real c,
                              rvec x[],rvec f[],rvec fshift[],
+                             gmx_bool VirCorr,matrix dxdf,
                              t_pbc *pbc,t_graph *g)
 {
     rvec xvi,xij,xik,xil,ra,rb,rja,rjb,rab,rm,rt;
@@ -1135,6 +1218,22 @@ static void spread_vsite4FDN(t_iatom ia[],real a,real b,real c,
         rvec_inc(fshift[sik],fk);
         rvec_inc(fshift[sil],fl);
     }
+
+    if (VirCorr)
+    {
+        rvec xiv;
+        int  i,j;
+
+        pbc_rvec_sub(pbc,x[av],x[ai],xiv);
+
+        for(i=0; i<DIM; i++)
+        {
+            for(j=0; j<DIM; j++)
+            {
+                dxdf[i][j] += -xiv[i]*fv[j] + xij[i]*fj[j] + xik[i]*fk[j] + xil[i]*fl[j];
+            }
+        }
+    }
     
     /* Total: 207 flops (Yuck!) */
 }
@@ -1179,10 +1278,11 @@ static int spread_vsiten(t_iatom ia[],t_iparams ip[],
 
 
 void spread_vsite_f(FILE *log,gmx_vsite_t *vsite,
-                   rvec x[],rvec f[],rvec *fshift,
-                   t_nrnb *nrnb,t_idef *idef,
-                   int ePBC,gmx_bool bMolPBC,t_graph *g,matrix box,
-                   t_commrec *cr)
+                    rvec x[],rvec f[],rvec *fshift,
+                    gmx_bool VirCorr,matrix vir,
+                    t_nrnb *nrnb,t_idef *idef,
+                    int ePBC,gmx_bool bMolPBC,t_graph *g,matrix box,
+                    t_commrec *cr)
 {
   real      a1,b1,c1;
   int       i,inc,m,nra,nr,tp,ftype;
@@ -1191,6 +1291,7 @@ void spread_vsite_f(FILE *log,gmx_vsite_t *vsite,
   t_iparams *ip;
   t_pbc     pbc,*pbc_null,*pbc_null2;
   int       *vsite_pbc;
+  matrix    dxdf;
 
   /* We only need to do pbc when we have inter-cg vsites */
   if ((DOMAINDECOMP(cr) || bMolPBC) && vsite->n_intercg_vsite) {
@@ -1210,6 +1311,11 @@ void spread_vsite_f(FILE *log,gmx_vsite_t *vsite,
   {
     pd_clear_nonlocal_constructs(vsite->vsitecomm,f);
   }
+
+    if (vir != NULL)
+    {
+        clear_mat(dxdf);
+    }
        
   ip     = idef->iparams;
 
@@ -1266,30 +1372,30 @@ void spread_vsite_f(FILE *log,gmx_vsite_t *vsite,
          break;
        case F_VSITE3FD:
          b1 = ip[tp].vsite.b;
-         spread_vsite3FD(ia,a1,b1,x,f,fshift,pbc_null2,g);
+         spread_vsite3FD(ia,a1,b1,x,f,fshift,VirCorr,dxdf,pbc_null2,g);
          nd3FD++;
          break;
        case F_VSITE3FAD:
          b1 = ip[tp].vsite.b;
-         spread_vsite3FAD(ia,a1,b1,x,f,fshift,pbc_null2,g);
+         spread_vsite3FAD(ia,a1,b1,x,f,fshift,VirCorr,dxdf,pbc_null2,g);
          nd3FAD++;
          break;
        case F_VSITE3OUT:
          b1 = ip[tp].vsite.b;
          c1 = ip[tp].vsite.c;
-         spread_vsite3OUT(ia,a1,b1,c1,x,f,fshift,pbc_null2,g);
+         spread_vsite3OUT(ia,a1,b1,c1,x,f,fshift,VirCorr,dxdf,pbc_null2,g);
          nd3OUT++;
          break;
        case F_VSITE4FD:
          b1 = ip[tp].vsite.b;
          c1 = ip[tp].vsite.c;
-         spread_vsite4FD(ia,a1,b1,c1,x,f,fshift,pbc_null2,g);
+         spread_vsite4FD(ia,a1,b1,c1,x,f,fshift,VirCorr,dxdf,pbc_null2,g);
          nd4FD++;
          break;
        case F_VSITE4FDN:
          b1 = ip[tp].vsite.b;
          c1 = ip[tp].vsite.c;
-         spread_vsite4FDN(ia,a1,b1,c1,x,f,fshift,pbc_null2,g);
+         spread_vsite4FDN(ia,a1,b1,c1,x,f,fshift,VirCorr,dxdf,pbc_null2,g);
          nd4FDN++;
          break;
        case F_VSITEN:
@@ -1308,7 +1414,20 @@ void spread_vsite_f(FILE *log,gmx_vsite_t *vsite,
       }
     }
   }
-       
+
+    if (VirCorr)
+    {
+        int i,j;
+
+        for(i=0; i<DIM; i++)
+        {
+            for(j=0; j<DIM; j++)
+            {
+                vir[i][j] += -0.5*dxdf[i][j];
+            }
+        }
+    }
+
   inc_nrnb(nrnb,eNR_VSITE2,   nd2     );
   inc_nrnb(nrnb,eNR_VSITE3,   nd3     );
   inc_nrnb(nrnb,eNR_VSITE3FD, nd3FD   );

diff --git a/src/tools/gmx_bar.c b/src/tools/gmx_bar.c
index e5fdc0335bfa895893bec8ec42962c572dfe8146..a82661235d87d11557b941c2f5fc91ccef55d609 100644 (file)
--- a/src/tools/gmx_bar.c
+++ b/src/tools/gmx_bar.c
@@ -2030,7 +2030,7 @@ static void read_bar_xvg(char *fn, real *temp, lambda_t *lambda_head)
         gmx_fatal(FARGS,"File '%s' contains fewer than two columns", fn);
     }
 
-    if ( ( *temp != barsim->temp) && (*temp > 0) )
+    if ( !gmx_within_tol(*temp,barsim->temp,GMX_FLOAT_EPS) && (*temp > 0) )
     {
         gmx_fatal(FARGS,"Temperature in file %s different from earlier files or setting\n", fn);
     }