From a2367f3cbab136cc0f521500032d7fd652590625 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Mon, 29 Jul 2013 17:32:29 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@10466
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/USER-COLVARS/fix_colvars.cpp | 12 +++++++++++-
 src/dump_image.cpp               |  4 ++++
 src/fix_shake.cpp                | 27 +++++++++++++++++++--------
 src/neigh_derive.cpp             |  4 ++++
 src/respa.cpp                    |  3 ++-
 5 files changed, 40 insertions(+), 10 deletions(-)

diff --git a/src/USER-COLVARS/fix_colvars.cpp b/src/USER-COLVARS/fix_colvars.cpp
index f93c6c1156..3079b0e8ad 100644
--- a/src/USER-COLVARS/fix_colvars.cpp
+++ b/src/USER-COLVARS/fix_colvars.cpp
@@ -35,6 +35,16 @@
 
 #include "colvarproxy_lammps.h"
 
+static const char colvars_pub[] =
+  "@Article{fiorin13,\n"
+  " author =  {G.~Fiorin and M.{\\,}L.~Klein and J.~H{\\'e}nin},\n"
+  " title =   {Using collective variables to drive molecular"
+  " dynamics simulations},\n"
+  " journal = {Mol.~Phys.},\n"
+  " year =    2013,\n"
+  " note =    {doi: 10.1080/00268976.2013.813594}\n"
+  "}\n\n";
+
 /* re-usable integer hash table code with static linkage. */
 
 /** hash table top level data structure */
@@ -360,7 +370,7 @@ FixColvars::FixColvars(LAMMPS *lmp, int narg, char **arg) :
   /* storage required to communicate a single coordinate or force. */
   size_one = sizeof(struct commdata);
 
-  citeme->add(CiteMe::FIORIN_2013);
+  citeme->add(colvars_pub);
 }
 
 /*********************************
diff --git a/src/dump_image.cpp b/src/dump_image.cpp
index 307eb927e3..2f302354b8 100644
--- a/src/dump_image.cpp
+++ b/src/dump_image.cpp
@@ -47,6 +47,10 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) :
 {
   if (binary || multiproc) error->all(FLERR,"Invalid dump image filename");
 
+  // force binary flag on to avoid corrupted output on Windows
+
+  binary = 1;
+
   // set filetype based on filename suffix
 
   int n = strlen(filename);
diff --git a/src/fix_shake.cpp b/src/fix_shake.cpp
index a66c8d554a..164614972a 100644
--- a/src/fix_shake.cpp
+++ b/src/fix_shake.cpp
@@ -413,9 +413,15 @@ void FixShake::setup(int vflag)
     dtv = step_respa[0];
     dtf_innerhalf = 0.5 * step_respa[0] * force->ftm2v;
     dtf_inner = dtf_innerhalf;
-    ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
-    post_force_respa(vflag,nlevels_respa-1,0);
-    ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
+
+    // apply correction to all rRESPA levels
+
+    for (int ilevel = 0; ilevel < nlevels_respa; ilevel++) {
+      ((Respa *) update->integrate)->copy_flevel_f(ilevel);
+      post_force_respa(vflag,ilevel,loop_respa[ilevel]-1);
+      ((Respa *) update->integrate)->copy_f_flevel(ilevel);
+    }
+
     dtf_inner = step_respa[0] * force->ftm2v;
   }
 }
@@ -543,10 +549,12 @@ void FixShake::post_force_respa(int vflag, int ilevel, int iloop)
 
   if (ilevel == nlevels_respa-1 && update->ntimestep == next_output) stats();
 
-  // enforce SHAKE constraints on every loop iteration of every rRESPA level
-  // except last loop iteration of inner levels
+  // might be OK to skip enforcing SHAKE constraings
+  // on last iteration of inner levels if pressure not requested
+  // however, leads to slightly different trajectories
 
-  if (ilevel < nlevels_respa-1 && iloop == loop_respa[ilevel]-1) return;
+  //if (ilevel < nlevels_respa-1 && iloop == loop_respa[ilevel]-1 && !vflag)
+  //  return;
 
   // xshake = unconstrained move with current v,f as function of level
   // communicate results if necessary
@@ -554,9 +562,12 @@ void FixShake::post_force_respa(int vflag, int ilevel, int iloop)
   unconstrained_update_respa(ilevel);
   if (nprocs > 1) comm->forward_comm_fix(this);
 
-  // virial setup, only need to compute on outermost level
+  // virial setup only needed on last iteration of innermost level
+  //   and if pressure is requested
+  // virial accumulation happens at outermost level
 
-  if (ilevel == nlevels_respa-1 && vflag) v_setup(vflag);
+  if (ilevel == 0 && iloop == loop_respa[ilevel] - 1 && vflag) v_setup(vflag);
+  if (iloop == loop_respa[ilevel]-1) evflag = 1;
   else evflag = 0;
 
   // loop over clusters to add constraint forces
diff --git a/src/neigh_derive.cpp b/src/neigh_derive.cpp
index dbf5fe258a..5de3750968 100644
--- a/src/neigh_derive.cpp
+++ b/src/neigh_derive.cpp
@@ -439,9 +439,13 @@ void Neighbor::skip_from_respa(NeighList *list)
       if (ipage_middle->status())
         error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
     }
+
+    inum++;
   }
 
   list->inum = inum;
+  listinner->inum = inum;
+  if (respamiddle) listmiddle->inum = inum;
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/respa.cpp b/src/respa.cpp
index 67783ab2d1..19a6f7d800 100644
--- a/src/respa.cpp
+++ b/src/respa.cpp
@@ -286,7 +286,8 @@ void Respa::init()
 
   ev_setup();
 
-  // detect if fix omp is present and will clear force arrays for us
+  // detect if fix omp is present and will clear force arrays
+
   int ifix = modify->find_fix("package_omp");
   if (ifix >= 0) external_force_clear = 1;