Merge "fixed initial COM motion removal with multiple groups" into release-4-5-patches

diff --git a/src/kernel/md.c b/src/kernel/md.c
index 67e51a7a0d415f0f09d363c67e57dc1ebebf68bf..c5dd33da726eb6538a5f982390a9e56998b96f24 100644 (file)
--- a/src/kernel/md.c
+++ b/src/kernel/md.c
@@ -469,6 +469,7 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
   
     /* I'm assuming we need global communication the first time! MRS */
     cglo_flags = (CGLO_TEMPERATURE | CGLO_GSTAT
+                  | ((ir->comm_mode != ecmNO) ? CGLO_STOPCM:0)
                   | (bVV ? CGLO_PRESSURE:0)
                   | (bVV ? CGLO_CONSTRAINT:0)
                   | (bRerunMD ? CGLO_RERUNMD:0)
@@ -490,7 +491,7 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
                         NULL,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
                         constr,NULL,FALSE,state->box,
                         top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
-                        cglo_flags &~ CGLO_PRESSURE);
+                        cglo_flags &~ (CGLO_STOPCM | CGLO_PRESSURE));
     }
     
     /* Calculate the initial half step temperature, and save the ekinh_old */
@@ -972,7 +973,6 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
         
         /* these CGLO_ options remain the same throughout the iteration */
         cglo_flags = ((bRerunMD ? CGLO_RERUNMD : 0) |
-                      (bStopCM ? CGLO_STOPCM : 0) |
                       (bGStat ? CGLO_GSTAT : 0)
             );
         
@@ -1133,6 +1133,7 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
                                 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
                                 cglo_flags 
                                 | CGLO_ENERGY 
+                                | (bStopCM ? CGLO_STOPCM : 0)
                                 | (bTemp ? CGLO_TEMPERATURE:0) 
                                 | (bPres ? CGLO_PRESSURE : 0) 
                                 | (bPres ? CGLO_CONSTRAINT : 0)
@@ -1574,6 +1575,7 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
                             top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
                             cglo_flags 
                             | (!EI_VV(ir->eI) ? CGLO_ENERGY : 0) 
+                            | (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0)
                             | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0) 
                             | (!EI_VV(ir->eI) || bRerunMD ? CGLO_PRESSURE : 0) 
                             | (bIterations && iterate.bIterate ? CGLO_ITERATE : 0) 

diff --git a/src/mdlib/md_support.c b/src/mdlib/md_support.c
index 7551a76d8c2b3213b574f02d7d7961962ffd8f0c..edad6659cf3c400788421332e600bd021d84f113 100644 (file)
--- a/src/mdlib/md_support.c
+++ b/src/mdlib/md_support.c
@@ -342,16 +342,16 @@ void compute_globals(FILE *fplog, gmx_global_stat_t gstat, t_commrec *cr, t_inpu
         }
         
         debug_gmx();
-        
-        /* Calculate center of mass velocity if necessary, also parallellized */
-        if (bStopCM && !bRerunMD && bEner) 
-        {
-            calc_vcm_grp(fplog,mdatoms->start,mdatoms->homenr,mdatoms,
-                         state->x,state->v,vcm);
-        }
     }
 
-    if (bTemp || bPres || bEner || bConstrain) 
+    /* Calculate center of mass velocity if necessary, also parallellized */
+    if (bStopCM)
+    {
+        calc_vcm_grp(fplog,mdatoms->start,mdatoms->homenr,mdatoms,
+                     state->x,state->v,vcm);
+    }
+
+    if (bTemp || bStopCM || bPres || bEner || bConstrain)
     {
         if (!bGStat)
         {
@@ -375,7 +375,7 @@ void compute_globals(FILE *fplog, gmx_global_stat_t gstat, t_commrec *cr, t_inpu
                 wallcycle_start(wcycle,ewcMoveE);
                 GMX_MPE_LOG(ev_global_stat_start);
                 global_stat(fplog,gstat,cr,enerd,force_vir,shake_vir,mu_tot,
-                            ir,ekind,constr,vcm,
+                            ir,ekind,constr,bStopCM ? vcm : NULL,
                             gs != NULL ? eglsNR : 0,gs_buf,
                             top_global,state,
                             *bSumEkinhOld,flags);
@@ -425,16 +425,17 @@ void compute_globals(FILE *fplog, gmx_global_stat_t gstat, t_commrec *cr, t_inpu
                      mdatoms->massT,mdatoms->tmass,ekind->ekin);
     }
     
-    if (bEner) {
-        /* Do center of mass motion removal */
-        if (bStopCM && !bRerunMD) /* is this correct?  Does it get called too often with this logic? */
-        {
-            check_cm_grp(fplog,vcm,ir,1);
-            do_stopcm_grp(fplog,mdatoms->start,mdatoms->homenr,mdatoms->cVCM,
-                          state->x,state->v,vcm);
-            inc_nrnb(nrnb,eNR_STOPCM,mdatoms->homenr);
-        }
+    /* Do center of mass motion removal */
+    if (bStopCM)
+    {
+        check_cm_grp(fplog,vcm,ir,1);
+        do_stopcm_grp(fplog,mdatoms->start,mdatoms->homenr,mdatoms->cVCM,
+                      state->x,state->v,vcm);
+        inc_nrnb(nrnb,eNR_STOPCM,mdatoms->homenr);
+    }
 
+    if (bEner)
+    {
         /* Calculate the amplitude of the cosine velocity profile */
         ekind->cosacc.vcos = ekind->cosacc.mvcos/mdatoms->tmass;
     }

diff --git a/src/mdlib/sim_util.c b/src/mdlib/sim_util.c
index 590aaaf5672ac4d7507a00e1ddffb8d57784ecdf..2d7fd82a742210b9b24fa90bb9382f07e89f6ed4 100644 (file)
--- a/src/mdlib/sim_util.c
+++ b/src/mdlib/sim_util.c
@@ -921,7 +921,6 @@ void do_constrain_first(FILE *fplog,gmx_constr_t constr,
 {
     int    i,m,start,end;
     gmx_large_int_t step;
-    double mass,tmass,vcm[4];
     real   dt=ir->delta_t;
     real   dvdlambda;
     rvec   *savex;
@@ -998,38 +997,6 @@ void do_constrain_first(FILE *fplog,gmx_constr_t constr,
         }
     }
     
-    for(m=0; (m<4); m++)
-        vcm[m] = 0;
-    for(i=start; i<end; i++) {
-        mass = md->massT[i];
-        for(m=0; m<DIM; m++) {
-            vcm[m] += state->v[i][m]*mass;
-        }
-        vcm[3] += mass;
-    }
-    
-    if (ir->nstcomm != 0 || debug) {
-        /* Compute the global sum of vcm */
-        if (debug)
-            fprintf(debug,"vcm: %8.3f  %8.3f  %8.3f,"
-                    " total mass = %12.5e\n",vcm[XX],vcm[YY],vcm[ZZ],vcm[3]);
-        if (PAR(cr))
-            gmx_sumd(4,vcm,cr);
-        tmass = vcm[3];
-        for(m=0; (m<DIM); m++)
-            vcm[m] /= tmass;
-        if (debug) 
-            fprintf(debug,"vcm: %8.3f  %8.3f  %8.3f,"
-                    " total mass = %12.5e\n",vcm[XX],vcm[YY],vcm[ZZ],tmass);
-        if (ir->nstcomm != 0) {
-            /* Now we have the velocity of center of mass, let's remove it */
-            for(i=start; (i<end); i++) {
-                for(m=0; (m<DIM); m++)
-                    state->v[i][m] -= vcm[m];
-            }
-
-        }
-    }
     sfree(savex);
 }
 

diff --git a/src/mdlib/stat.c b/src/mdlib/stat.c
index 185b15cb6b81d7c62436c9d34f39df985d347229..d84ab3f2b61a1fee3171b269700d07ee840edfb6 100644 (file)
--- a/src/mdlib/stat.c
+++ b/src/mdlib/stat.c
@@ -273,24 +273,25 @@ void global_stat(FILE *fplog,gmx_global_stat_t gs,
               }
           }
       }
+  }
 
-      if (vcm) 
+  if (vcm)
+  {
+      icm   = add_binr(rb,DIM*vcm->nr,vcm->group_p[0]);
+      where();
+      imass = add_binr(rb,vcm->nr,vcm->group_mass);
+      where();
+      if (vcm->mode == ecmANGULAR)
       {
-          icm   = add_binr(rb,DIM*vcm->nr,vcm->group_p[0]);
+          icj   = add_binr(rb,DIM*vcm->nr,vcm->group_j[0]);
           where();
-          imass = add_binr(rb,vcm->nr,vcm->group_mass);
+          icx   = add_binr(rb,DIM*vcm->nr,vcm->group_x[0]);
+          where();
+          ici   = add_binr(rb,DIM*DIM*vcm->nr,vcm->group_i[0][0]);
           where();
-          if (vcm->mode == ecmANGULAR) 
-          {
-              icj   = add_binr(rb,DIM*vcm->nr,vcm->group_j[0]);
-              where();
-              icx   = add_binr(rb,DIM*vcm->nr,vcm->group_x[0]);
-              where();
-              ici   = add_binr(rb,DIM*DIM*vcm->nr,vcm->group_i[0][0]);
-              where();
-          }
       }
   }
+
   if (DOMAINDECOMP(cr)) 
   {
       nb = cr->dd->nbonded_local;
@@ -373,23 +374,6 @@ void global_stat(FILE *fplog,gmx_global_stat_t gs,
                   extract_bind(rb,iepl,enerd->n_lambda,enerd->enerpart_lambda);
               }
           }
-          /* should this be here, or with ekin?*/
-          if (vcm) 
-          {
-              extract_binr(rb,icm,DIM*vcm->nr,vcm->group_p[0]);
-              where();
-              extract_binr(rb,imass,vcm->nr,vcm->group_mass);
-              where();
-              if (vcm->mode == ecmANGULAR) 
-              {
-                  extract_binr(rb,icj,DIM*vcm->nr,vcm->group_j[0]);
-                  where();
-                  extract_binr(rb,icx,DIM*vcm->nr,vcm->group_x[0]);
-                  where();
-                  extract_binr(rb,ici,DIM*DIM*vcm->nr,vcm->group_i[0][0]);
-                  where();
-              }
-          }
           if (DOMAINDECOMP(cr)) 
           {
               extract_bind(rb,inb,1,&nb);
@@ -406,6 +390,23 @@ void global_stat(FILE *fplog,gmx_global_stat_t gs,
       }
   }
 
+  if (vcm)
+  {
+      extract_binr(rb,icm,DIM*vcm->nr,vcm->group_p[0]);
+      where();
+      extract_binr(rb,imass,vcm->nr,vcm->group_mass);
+      where();
+      if (vcm->mode == ecmANGULAR)
+      {
+          extract_binr(rb,icj,DIM*vcm->nr,vcm->group_j[0]);
+          where();
+          extract_binr(rb,icx,DIM*vcm->nr,vcm->group_x[0]);
+          where();
+          extract_binr(rb,ici,DIM*DIM*vcm->nr,vcm->group_i[0][0]);
+          where();
+      }
+  }
+
   if (nsig > 0) 
   {
       extract_binr(rb,isig,nsig,sig);