git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@15019 f3b2605a-c512-4ea7-a41b-209d697bcdaa

src/RIGID/fix_ehex.cpp

@@ -0,0 +1,625 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author: Peter Wirnsberger (University of Cambridge)
+
+   This source file implements the asymmetric version of the enhanced heat
+   exchange (eHEX/a) algorithm. The paper is available for download on
+   arXiv: http://arxiv.org/pdf/1507.07081.pdf.
+
+   This file is based on fix_heat.cpp written by Paul Crozier (SNL)
+   which implements the heat exchange (HEX) algorithm.
+------------------------------------------------------------------------- */
+
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+#include "fix_ehex.h"
+#include "atom.h"
+#include "domain.h"
+#include "region.h"
+#include "group.h"
+#include "force.h"
+#include "update.h"
+#include "modify.h"
+#include "input.h"
+#include "variable.h"
+#include "memory.h"
+#include "error.h"
+#include "fix_shake.h"
+#include "neighbor.h"
+#include "comm.h"
+#include "timer.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+enum{CONSTANT,EQUAL,ATOM};
+
+/* ---------------------------------------------------------------------- */
+
+FixEHEX::FixEHEX(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
+{
+  MPI_Comm_rank(world, &me);
+ 
+  // check
+  if (narg < 4) error->all(FLERR,"Illegal fix ehex command: wrong number of parameters ");
+
+  scalar_flag = 1;
+  global_freq = 1;
+  extscalar = 0;
+
+  // apply fix every nevery timesteps
+
+  nevery = force->inumeric(FLERR,arg[3]);
+
+  if (nevery <= 0) error->all(FLERR,"Illegal fix ehex command");
+
+  // heat flux into the reservoir
+
+  heat_input = force->numeric(FLERR,arg[4]);
+
+  // optional args
+
+  iregion = -1;
+  idregion = NULL;
+
+  // NOTE: constraints are deactivated by default
+
+  constraints = 0;
+
+  // NOTE: cluster rescaling is deactivated by default
+
+  cluster = 0;
+
+  // NOTE: hex = 1 means that no coordinate correction is applied in which case eHEX reduces to HEX
+ 
+  hex = 0;
+
+  int iarg = 5;
+  while (iarg < narg) {
+
+    if (strcmp(arg[iarg],"region") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal fix ehex command: wrong number of parameters ");
+      iregion = domain->find_region(arg[iarg+1]);
+      if (iregion == -1)
+        error->all(FLERR,"Region ID for fix ehex does not exist");
+      int n = strlen(arg[iarg+1]) + 1;
+      idregion = new char[n];
+      strcpy(idregion,arg[iarg+1]);
+      iarg += 2;
+    } 
+
+    // apply constraints (shake/rattle) at the end of the timestep
+
+    else if (strcmp(arg[iarg], "constrain") == 0) {
+      constraints = 1;
+      iarg += 1;
+    }
+
+    // rescale only if the entire molecule is contained within the region
+
+    else if (strcmp(arg[iarg], "com") == 0) {
+      cluster = 1;
+      iarg += 1;
+    }
+    
+    // don't apply a coordinate correction if this keyword is specified
+
+    else if (strcmp(arg[iarg], "hex") == 0) {
+      hex  = 1;
+      iarg+= 1;
+    }
+    else
+      error->all(FLERR, "Illegal fix ehex keyword ");
+  } 
+
+  // check options
+
+  if (cluster && !constraints) 
+    error->all(FLERR, "You can only use the keyword 'com' together with the keyword 'constrain' ");
+
+  scale = 1.0;
+  scalingmask    = NULL;
+  grow_arrays(atom->nmax);
+  atom->add_callback(0);
+
+}
+
+
+/* ---------------------------------------------------------------------- */
+
+void FixEHEX::grow_arrays(int nmax) {
+  memory->grow(scalingmask, nmax,"ehex:scalingmask");
+}
+
+
+/* ---------------------------------------------------------------------- */
+
+FixEHEX::~FixEHEX()
+{
+  atom->delete_callback(id,0);
+  delete [] idregion;
+  memory->destroy(scalingmask);
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixEHEX::setmask()
+{
+  int mask = 0;
+  mask |= END_OF_STEP;
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixEHEX::init()
+{
+  // set index and check validity of region
+
+  if (iregion >= 0) {
+    iregion = domain->find_region(idregion);
+    if (iregion == -1)
+      error->all(FLERR,"Region ID for fix ehex does not exist");
+  }
+
+  // cannot have 0 atoms in group
+
+  if (group->count(igroup) == 0)
+    error->all(FLERR,"Fix ehex group has no atoms");
+
+  fshake = NULL;
+  if (constraints) { 
+
+    // check if constraining algorithm is used (FixRattle inherits from FixShake)
+
+    int cnt_shake = 0;
+    int id_shake;
+    for (int i = 0; i < modify->nfix; i++) { 
+      if (strcmp("rattle", modify->fix[i]->style) == 0 ||
+          strcmp("shake", modify->fix[i]->style) == 0) {
+        cnt_shake++;
+        id_shake = i;
+      }
+    } 
+ 
+    if (cnt_shake > 1)         
+      error->all(FLERR,"Multiple instances of fix shake/rattle detected (not supported yet)");
+    else if (cnt_shake == 1)   {
+     fshake = ((FixShake*) modify->fix[id_shake]);
+    }
+    else if (cnt_shake == 0)                  
+      error->all(FLERR, "Fix ehex was configured with keyword constrain, but shake/rattle was not defined");
+  }
+}
+
+
+
+/* ---------------------------------------------------------------------- */
+
+
+void FixEHEX::end_of_step() {
+  // store local pointers
+
+  x       = atom->x;
+  f       = atom->f;
+  v       = atom->v;
+  mass    = atom->mass;
+  rmass   = atom->rmass;
+  type    = atom->type;
+  nlocal  = atom->nlocal;
+
+  // determine which sites are to be rescaled
+
+  update_scalingmask();
+
+  // rescale velocities 
+
+  rescale();
+
+  // if required use shake/rattle to correct coordinates and velocities
+
+  if (constraints && fshake) 
+    fshake->shake_end_of_step(0);
+}
+
+
+
+/* ----------------------------------------------------------------------
+   Iterate over all atoms, rescale the velocities and apply coordinate 
+   corrections.
+------------------------------------------------------------------------- */
+
+void FixEHEX::rescale() {
+  double heat,ke,massone;
+  double Kr, Ke, escale;
+  double vsub[3],vcm[3], sfr[3];
+  double dvcmdt[3];
+  double mi;
+  double dt;
+  double F, mr, Fo2Kr, epsr_ik, sfvr, eta_ik;
+
+  dt = update->dt;
+
+  // calculate centre of mass properties
+   
+  com_properties(vcm, sfr, &sfvr, &Ke, &Kr, &masstotal); 
+
+  // heat flux into the reservoir
+
+  F     = heat_input * force->ftm2v * nevery;
+  
+  // total mass
+
+  mr    = masstotal;
+  
+  Fo2Kr = F / (2.*Kr);
+
+  // energy scaling factor
+
+  escale = 1. + (F*dt)/Kr;
+
+  // safety check for kinetic energy
+
+  if (escale < 0.0) error->all(FLERR,"Fix ehex kinetic energy went negative");
+
+  scale = sqrt(escale);
+  vsub[0] = (scale-1.0) * vcm[0];
+  vsub[1] = (scale-1.0) * vcm[1];
+  vsub[2] = (scale-1.0) * vcm[2];
+
+  for (int i = 0; i < nlocal; i++){
+
+    if (scalingmask[i]) {
+
+      mi = (rmass) ? rmass[i] :  mass[type[i]];
+
+      for (int k=0; k<3; k++) {
+
+        // apply coordinate correction unless running in hex mode
+
+        if (!hex) {
+
+            // epsr_ik implements Eq. (20) in the paper
+
+            eta_ik    = mi * F/(2.*Kr) * (v[i][k] - vcm[k]);          
+            epsr_ik   = eta_ik / (mi*Kr) * (F/48. + sfvr/6.*force->ftm2v) - F/(12.*Kr) * (f[i][k]/mi - sfr[k]/mr)*force->ftm2v;
+
+            x[i][k]  -= dt*dt*dt * epsr_ik;
+        }
+
+        // rescale the velocity
+
+        v[i][k]   = scale*v[i][k] - vsub[k];
+      }
+    }
+  }
+}
+
+
+/* ---------------------------------------------------------------------- */
+
+double FixEHEX::compute_scalar()
+{
+  return scale;
+}
+
+
+/* ----------------------------------------------------------------------
+   memory usage of local atom-based arrays
+------------------------------------------------------------------------- */
+
+double FixEHEX::memory_usage()
+{
+  double bytes = 0.0;
+  bytes += atom->nmax * sizeof(double);
+  return bytes;
+}
+
+
+/* ----------------------------------------------------------------------
+   Update the array scalingmask depending on which individual atoms 
+   will be rescaled or not. 
+------------------------------------------------------------------------- */
+
+void FixEHEX::update_scalingmask() {
+  int m;
+  int lid;
+  bool stat;
+  int nsites;
+
+  // prematch region
+
+  Region *region = NULL;
+  if (iregion >= 0) {
+    region = domain->regions[iregion];
+    region->prematch();
+  }
+
+  // only rescale molecules whose center of mass if fully contained in the region
+
+  if (cluster) {
+ 
+    // loop over all clusters
+
+    for (int i=0; i < fshake->nlist; i++) {
+
+      // cluster id 
+
+      m    = fshake->list[i]; 
+
+      // check if the centre of mass of the cluster is inside the region
+      // if region == NULL, just check the group information of all sites
+
+      if      (fshake->shake_flag[m] == 1)      nsites = 3;
+      else if (fshake->shake_flag[m] == 2)      nsites = 2;
+      else if (fshake->shake_flag[m] == 3)      nsites = 3;
+      else if (fshake->shake_flag[m] == 4)      nsites = 4;
+      else                                      nsites = 0;
+
+      if (nsites == 0) {
+        error->all(FLERR,"Internal error: shake_flag[m] has to be between 1 and 4 for m in nlist");  
+      }
+
+      stat = check_cluster(fshake->shake_atom[m], nsites, region);
+
+      for (int l=0; l < nsites; l++)  {
+        lid = atom->map(fshake->shake_atom[m][l]);
+        scalingmask[lid] = stat;
+      }      
+    }
+
+    // check atoms that do not belong to any cluster
+
+    for (int i=0; i<atom->nlocal; i++)  {
+      if (fshake->shake_flag[i] == 0) 
+        scalingmask[i] = rescale_atom(i,region);
+    }
+
+  } 
+
+  // no clusters, just individual sites (e.g. monatomic system or flexible molecules)
+
+  else {
+    for (int i=0; i<atom->nlocal; i++) 
+      scalingmask[i] =  rescale_atom(i,region);
+  }
+
+}
+
+
+/* ----------------------------------------------------------------------
+   Check if the centre of mass of the cluster to be constrained is
+   inside the region.
+------------------------------------------------------------------------- */
+
+bool FixEHEX::check_cluster(int *shake_atom, int n, Region * region) {
+
+  // IMPORTANT NOTE: If any site of the cluster belongs to a group
+  //                 which should not be rescaled than all of the sites
+  //                 will be ignored!
+
+  double **x     = atom->x;
+  double * rmass = atom->rmass;
+  double * mass  = atom->mass; 
+  int    * type  = atom->type;
+  int    * mask  = atom->mask;
+  double   xcom[3], xtemp[3];
+  double   mcluster, mi;
+  bool     stat;
+  int      lid[4];
+
+  // accumulate mass and centre of mass position
+
+  stat      = true;
+  xcom[0]   = 0.;
+  xcom[1]   = 0.;
+  xcom[2]   = 0.;
+  mcluster  = 0;
+
+  for (int i = 0; i < n; i++) {
+
+    // get local id
+
+    lid[i] = atom->map(shake_atom[i]);
+    
+    // check if all sites of the cluster belong to the correct group
+
+    stat = stat && (mask[lid[i]] & groupbit);
+
+    if (region && stat)  {
+
+      // check if reduced mass is used
+
+      mi        = (rmass) ? rmass[lid[i]] : mass[type[lid[i]]];
+      mcluster += mi;
+
+      // accumulate centre of mass 
+      // NOTE: you can either use unwrapped coordinates or take site x[lid[0]] as reference,
+      //       i.e. reconstruct the molecule around this site and calculate the com.
+
+      for (int k=0; k<3; k++) 
+        xtemp[k] = x[lid[i]][k] - x[lid[0]][k];
+ 
+      // take into account pbc
+      
+      domain->minimum_image(xtemp);
+
+      for (int k=0; k<3; k++)
+        xcom[k] += mi * (x[lid[0]][k] + xtemp[k]) ;
+    }   
+  } 
+
+  // check if centre of mass is inside the region (if specified)
+
+  if (region && stat) {
+
+    // check mass
+
+    if (mcluster < 1.e-14) {
+      error->all(FLERR, "Fix ehex shake cluster has almost zero mass.");
+    }
+
+    // divide by total mass
+
+    for (int k=0; k<3; k++) 
+      xcom[k] = xcom[k]/mcluster;      
+ 
+    // apply periodic boundary conditions (centre of mass could be outside the box)
+    // and check if molecule is inside the region
+
+    domain->remap(xcom);
+    stat = stat && region->match(xcom[0], xcom[1], xcom[2]);
+  }
+  
+  return stat; 
+}
+
+
+/* ----------------------------------------------------------------------
+   Check if atom i has the correct group and is inside the region.
+------------------------------------------------------------------------- */
+
+bool FixEHEX::rescale_atom(int i, Region*region) {
+  bool stat;
+  double x_r[3];
+
+  // check mask and group
+
+  stat = (atom->mask[i] & groupbit);
+
+  if (region) {
+
+    x_r[0] = atom->x[i][0];
+    x_r[1] = atom->x[i][1];
+    x_r[2] = atom->x[i][2];
+
+    // apply periodic boundary conditions
+
+    domain->remap(x_r);
+
+    // check if the atom is in the group/region
+
+    stat = stat && region->match(x_r[0],x_r[1],x_r[2]);
+  }
+
+  return stat;
+}
+
+/* ----------------------------------------------------------------------
+   Calculate global properties of the atoms inside the reservoir.
+   (e.g. com velocity, kinetic energy, total mass,...)
+------------------------------------------------------------------------- */
+
+void FixEHEX::com_properties(double * vr, double * sfr, double *sfvr, double *K, double *Kr, double *mr) {
+   double ** f  = atom->f;
+   double ** v  = atom->v;
+   int nlocal   = atom->nlocal;
+   double *rmass= atom->rmass;
+   double *mass = atom->mass;
+   int    *type = atom->type;
+   double l_vr[3];
+   double l_mr;
+   double l_sfr[3];
+   double l_sfvr;
+   double l_K;
+   double mi;
+   double l_buf[9];
+   double buf[9];
+ 
+   // calculate partial sums
+
+   l_vr[0]  = l_vr[1]  = l_vr[2] = 0;
+   l_sfr[0] = l_sfr[1] = l_sfr[2] = 0;
+   l_sfvr   = 0;
+   l_mr     = 0;
+   l_K      = 0;
+
+   for (int i = 0; i < nlocal; i++) {
+     if (scalingmask[i]) {        
+       
+        // check if reduced mass is used
+
+        mi    = (rmass) ? rmass[i] : mass[type[i]];
+
+        // accumulate total mass
+
+        l_mr += mi;
+  
+        // accumulate kinetic energy
+
+        l_K  += mi/2. * (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]);
+
+        // sum_j f_j * v_j
+
+        l_sfvr  += f[i][0]*v[i][0] + f[i][1]*v[i][1] + f[i][2]*v[i][2];
+
+        // accumulate com velocity and sum of forces
+
+        for (int k=0; k<3; k++) {
+          l_vr[k] += mi * v[i][k];
+          l_sfr[k]+= f[i][k];
+        }
+     }
+   }
+
+   // reduce sums
+
+   l_buf[0] = l_vr[0];
+   l_buf[1] = l_vr[1];
+   l_buf[2] = l_vr[2];
+   l_buf[3] = l_K;
+   l_buf[4] = l_mr;
+   l_buf[5] = l_sfr[0];
+   l_buf[6] = l_sfr[1];
+   l_buf[7] = l_sfr[2];
+   l_buf[8] = l_sfvr;
+
+   MPI_Allreduce(l_buf, buf, 9, MPI_DOUBLE, MPI_SUM, world);
+
+   // total mass of region
+
+   *mr = buf[4];
+
+   if (*mr < 1.e-14) {
+      error->all(FLERR, "Fix ehex error mass of region is close to zero");
+   }
+
+   // total kinetic energy of region 
+
+   *K  = buf[3];
+
+   // centre of mass velocity of region
+
+   vr[0] = buf[0]/(*mr);
+   vr[1] = buf[1]/(*mr);
+   vr[2] = buf[2]/(*mr);
+
+   // sum of forces
+
+   sfr[0] = buf[5];
+   sfr[1] = buf[6];
+   sfr[2] = buf[7];
+
+   // calculate non-translational kinetic energy
+
+   *Kr = *K - 0.5* (*mr) * (vr[0]*vr[0]+vr[1]*vr[1]+vr[2]*vr[2]);
+
+   // calculate sum_j f_j * (v_j - v_r) = sum_j f_j * v_j  - v_r * sum_j f_j
+ 
+   *sfvr =  buf[8] - (vr[0]*sfr[0] + vr[1]*sfr[1] + vr[2]*sfr[2]);
+}
+

src/RIGID/fix_ehex.h

@@ -0,0 +1,130 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+
+#ifdef FIX_CLASS
+
+FixStyle(ehex,FixEHEX)
+
+#else
+
+#ifndef LMP_FIX_EHEX_H
+#define LMP_FIX_EHEX_H
+
+#include "fix.h"
+#include "fix_shake.h"
+#include "region.h"
+#define EHEX_DEBUG 0
+
+namespace LAMMPS_NS {
+
+class FixEHEX : public Fix {
+
+ public:
+  FixEHEX(class LAMMPS *, int, char **);
+  ~FixEHEX();
+  int  setmask();
+  void init();
+  void end_of_step();
+  void   rescale();
+  double compute_scalar();
+  double memory_usage();
+  void update_scalingmask();
+  void com_properties(double *, double *, double *, double*, double *, double*);
+  bool rescale_atom(int i, Region*region);
+  virtual void grow_arrays(int nmax);
+  bool check_cluster(int *shake_atom, int n, Region * region);
+
+ private:
+  int iregion;
+  double heat_input;
+  double masstotal;
+  double scale;
+  char *idregion;
+  int me;
+ 
+  double ** x;              // coordinates
+  double ** f;              // forces
+  double ** v;              // velocities
+  double *  mass;           // masses
+  double *  rmass;          // reduced masses
+  int    *  type;           // atom types
+  int   nlocal;             // number of local atoms
+  FixShake * fshake;        // pointer to fix_shake/fix_rattle
+  int constraints;          // constraints (0/1)
+  int cluster;              // rescaling entire clusters (0/1)
+  int hex;                  // HEX mode (0/1)
+  bool * scalingmask;       // scalingmask[i] determines whether 
+                            // the velocity of atom i is to be rescaled
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal fix ehex command: wrong number of parameters 
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Illegal fix ehex command: integer value expected
+
+Self-explanatory. Check the value for nevery.
+
+E: Region ID for fix ehex does not exist
+
+Self-explanatory.
+
+E: You can only use the keyword 'com' together with the keyword 'constrain' .
+
+Self-explanatory.
+
+E: Illegal fix ehex keyword
+
+Self-explanatory.
+
+E: Fix ehex group has no atoms
+
+Self-explanatory.
+
+E: Multiple instances of fix shake/rattle detected (not supported yet)
+
+You can only have one instance of fix rattle/shake at the moment.
+
+E: Fix ehex was configured with keyword constrain, but shake/rattle was not defined
+
+The option constrain requires either fix shake or fix rattle which is missing in the input script.
+
+E: Fix heat kinetic energy went negative
+
+This will cause the velocity rescaling about to be performed by fix
+heat to be invalid.
+
+E: Fix heat kinetic energy of an atom went negative
+
+This will cause the velocity rescaling about to be performed by fix
+heat to be invalid.
+
+E: Internal error: shake_flag[m] has to be between 1 and 4 for m in nlist
+
+Contact developers.
+
+E: Fix ehex error mass of region is close to zero
+
+Check your configuration.
+
+*/

src/RIGID/fix_rattle.cpp

@@ -44,13 +44,22 @@ using namespace FixConst;
 using namespace MathConst;
 using namespace MathExtra;
 
-// set RATTLE_DEBUG = 1 to check constraints at end of timestep
+// set RATTLE_DEBUG  1 to check constraints at end of timestep
 
 #define RATTLE_DEBUG 0
-#define RATTLE_TEST_VEL true
-#define RATTLE_TEST_POS true
 
-enum{V,VP,XSHAKE,X};
+// set RATTLE_RAISE_ERROR 1 if you want this fix to raise
+//     an error if the constraints cannot be satisfied
+
+#define RATTLE_RAISE_ERROR 0
+
+// You can enable velocity and coordinate checks separately
+// by setting RATTLE_TEST_VEL/POS true
+
+#define RATTLE_TEST_VEL false
+#define RATTLE_TEST_POS false
+
+enum{V,VP,XSHAKE};
 
 /* ---------------------------------------------------------------------- */
 
@@ -59,10 +68,6 @@ FixRattle::FixRattle(LAMMPS *lmp, int narg, char **arg) :
 {
   rattle = 1;
 
-  // define timestep for velocity integration
-
-  dtfv = 0.5 * update->dt * force->ftm2v;
-
   // allocate memory for unconstrained velocity update
 
   vp = NULL;
@@ -74,6 +79,9 @@ FixRattle::FixRattle(LAMMPS *lmp, int narg, char **arg) :
 
   comm_mode = XSHAKE;
   vflag_post_force = 0;
+
+  verr_max = 0;  
+  derr_max = 0;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -81,6 +89,25 @@ FixRattle::FixRattle(LAMMPS *lmp, int narg, char **arg) :
 FixRattle::~FixRattle()
 { 
   memory->destroy(vp);
+
+
+  if (RATTLE_DEBUG) {
+
+    // communicate maximum distance error
+
+    double global_derr_max, global_verr_max;
+    int npid;
+
+    MPI_Reduce(&derr_max, &global_derr_max, 1 , MPI_DOUBLE, MPI_MAX, 0, world);
+    MPI_Reduce(&verr_max, &global_verr_max, 1 , MPI_DOUBLE, MPI_MAX, 0, world);
+
+    MPI_Comm_rank (world, &npid); // Find out process rank
+
+    if (npid == 0 && screen) {
+      fprintf(screen, "RATTLE: Maximum overall relative position error ( (r_ij-d_ij)/d_ij ): %.10g\n", global_derr_max);
+      fprintf(screen, "RATTLE: Maximum overall absolute velocity error (r_ij * v_ij): %.10g\n", global_verr_max);
+    }
+  }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -118,7 +145,7 @@ void FixRattle::init() {
 
   if (flag && comm->me == 0)
     error->warning(FLERR,
-                   "Fix rattle should come after all other integration fixes");
+                   "Fix rattle should come after all other integration fixes ");
 }
 
 /* ----------------------------------------------------------------------
@@ -213,30 +240,36 @@ void FixRattle::final_integrate_respa(int ilevel, int iloop)
 void FixRattle::vrattle3angle(int m)
 {
   tagint i0,i1,i2;
+  int nlist,list[3];
   double c[3], l[3], a[3][3], r01[3], imass[3],
          r02[3], r12[3], vp01[3], vp02[3], vp12[3];
 
   // local atom IDs and constraint distances
+
   i0 = atom->map(shake_atom[m][0]);
   i1 = atom->map(shake_atom[m][1]);
   i2 = atom->map(shake_atom[m][2]);
 
   // r01,r02,r12 = distance vec between atoms
+
   MathExtra::sub3(x[i1],x[i0],r01);
   MathExtra::sub3(x[i2],x[i0],r02);
   MathExtra::sub3(x[i2],x[i1],r12);
 
   // take into account periodicity
+
   domain->minimum_image(r01);
   domain->minimum_image(r02);
   domain->minimum_image(r12);
 
   // v01,v02,v12 = velocity differences 
+
   MathExtra::sub3(vp[i1],vp[i0],vp01);
   MathExtra::sub3(vp[i2],vp[i0],vp02);
   MathExtra::sub3(vp[i2],vp[i1],vp12);
 
   // matrix coeffs and rhs for lamda equations
+
   if (rmass) {
     imass[0] = 1.0/rmass[i0];
     imass[1] = 1.0/rmass[i1];
@@ -248,6 +281,7 @@ void FixRattle::vrattle3angle(int m)
   }
 
   // setup matrix
+
   a[0][0]   =   (imass[1] + imass[0])   * MathExtra::dot3(r01,r01);
   a[0][1]   =   (imass[0]           )   * MathExtra::dot3(r01,r02);
   a[0][2]   =   (-imass[1]          )   * MathExtra::dot3(r01,r12);
@@ -259,14 +293,17 @@ void FixRattle::vrattle3angle(int m)
   a[2][2]   =   (imass[2] + imass[1])   * MathExtra::dot3(r12,r12);
 
   // sestup RHS
+
   c[0]  = -MathExtra::dot3(vp01,r01);
   c[1]  = -MathExtra::dot3(vp02,r02);
   c[2]  = -MathExtra::dot3(vp12,r12);
 
   // calculate the inverse matrix exactly
+
   solve3x3exactly(a,c,l);
 
   // add corrections to the velocities if processor owns atom
+
   if (i0 < nlocal) {
     for (int k=0; k<3; k++)  
       v[i0][k]  -=  imass[0]*  (  l[0] * r01[k] + l[1] * r02[k] );
@@ -286,20 +323,25 @@ void FixRattle::vrattle3angle(int m)
 void FixRattle::vrattle2(int m)
 {
   tagint    i0, i1;
+  int       nlist,list[2];
   double    imass[2], r01[3], vp01[3];
 
   // local atom IDs and constraint distances
+
   i0 = atom->map(shake_atom[m][0]);
   i1 = atom->map(shake_atom[m][1]);
 
   // r01 = distance vec between atoms, with PBC
+
   MathExtra::sub3(x[i1],x[i0],r01);
   domain->minimum_image(r01);
 
   // v01 = distance vectors for velocities
+
   MathExtra::sub3(vp[i1],vp[i0],vp01);
 
   // matrix coeffs and rhs for lamda equations
+
   if (rmass) {
     imass[0] = 1.0/rmass[i0];
     imass[1] = 1.0/rmass[i1];
@@ -309,10 +351,12 @@ void FixRattle::vrattle2(int m)
   }
 
   // Lagrange multiplier: exact solution
+
   double l01 = - MathExtra::dot3(r01,vp01) / 
     (MathExtra::dot3(r01,r01) * (imass[0] + imass[1]));
 
   // add corrections to the velocities if the process owns this atom
+
   if (i0 < nlocal) {
     for (int k=0; k<3; k++)  
       v[i0][k] -= imass[0] * l01 * r01[k];
@@ -328,15 +372,18 @@ void FixRattle::vrattle2(int m)
 void FixRattle::vrattle3(int m)
 {
   tagint    i0,i1,i2;
+  int       nlist,list[3];
   double    imass[3], r01[3], r02[3], vp01[3], vp02[3],
             a[2][2],c[2],l[2];
 
   // local atom IDs and constraint distances
+
   i0 = atom->map(shake_atom[m][0]);
   i1 = atom->map(shake_atom[m][1]);
   i2 = atom->map(shake_atom[m][2]);
 
   // r01,r02 = distance vec between atoms, with PBC
+
   MathExtra::sub3(x[i1],x[i0],r01);
   MathExtra::sub3(x[i2],x[i0],r02);
  
@@ -344,6 +391,7 @@ void FixRattle::vrattle3(int m)
   domain->minimum_image(r02);
 
   // vp01,vp02 =  distance vectors between velocities
+
   MathExtra::sub3(vp[i1],vp[i0],vp01);
   MathExtra::sub3(vp[i2],vp[i0],vp02);
 
@@ -358,19 +406,23 @@ void FixRattle::vrattle3(int m)
   }
 
   // setup matrix
+
   a[0][0]   =   (imass[1] + imass[0])   * MathExtra::dot3(r01,r01);
   a[0][1]   =   (imass[0]           )   * MathExtra::dot3(r01,r02);
   a[1][0]   =   a[0][1];
   a[1][1]   =   (imass[0] + imass[2])   * MathExtra::dot3(r02,r02);
 
   // setup RHS
+
   c[0]  = - MathExtra::dot3(vp01,r01);
   c[1]  = - MathExtra::dot3(vp02,r02);
 
   // calculate the inverse 2x2 matrix exactly
+
   solve2x2exactly(a,c,l);
 
   // add corrections to the velocities if the process owns this atom
+
   if (i0 < nlocal) {
     for (int k=0; k<3; k++)  
       v[i0][k] -= imass[0] * (  l[0] * r01[k] + l[1] * r02[k] );
@@ -390,16 +442,19 @@ void FixRattle::vrattle3(int m)
 void FixRattle::vrattle4(int m)
 {
   tagint    i0,i1,i2,i3;
+  int       nlist,list[4];
   double    imass[4], c[3], l[3], a[3][3],
             r01[3], r02[3], r03[3], vp01[3], vp02[3], vp03[3];
 
   // local atom IDs and constraint distances
+
   i0 = atom->map(shake_atom[m][0]);
   i1 = atom->map(shake_atom[m][1]);
   i2 = atom->map(shake_atom[m][2]);
   i3 = atom->map(shake_atom[m][3]);
 
   // r01,r02,r12 = distance vec between atoms, with PBC
+
   MathExtra::sub3(x[i1],x[i0],r01);
   MathExtra::sub3(x[i2],x[i0],r02);
   MathExtra::sub3(x[i3],x[i0],r03);
@@ -409,11 +464,13 @@ void FixRattle::vrattle4(int m)
   domain->minimum_image(r03);
 
   // vp01,vp02,vp03 = distance vectors between velocities
+
   MathExtra::sub3(vp[i1],vp[i0],vp01);
   MathExtra::sub3(vp[i2],vp[i0],vp02);
   MathExtra::sub3(vp[i3],vp[i0],vp03);
 
   // matrix coeffs and rhs for lamda equations
+
   if (rmass) {
     imass[0] = 1.0/rmass[i0];
     imass[1] = 1.0/rmass[i1];
@@ -427,6 +484,7 @@ void FixRattle::vrattle4(int m)
   }
 
   // setup matrix
+
   a[0][0]   =   (imass[0] + imass[1])   * MathExtra::dot3(r01,r01);
   a[0][1]   =   (imass[0]           )   * MathExtra::dot3(r01,r02);
   a[0][2]   =   (imass[0]           )   * MathExtra::dot3(r01,r03);
@@ -438,14 +496,17 @@ void FixRattle::vrattle4(int m)
   a[2][2]   =   (imass[0] + imass[3])   * MathExtra::dot3(r03,r03);
 
   // setup RHS
+
   c[0]  = - MathExtra::dot3(vp01,r01);
   c[1]  = - MathExtra::dot3(vp02,r02);
   c[2]  = - MathExtra::dot3(vp03,r03);
 
   // calculate the inverse 3x3 matrix exactly
+
   solve3x3exactly(a,c,l);
 
   // add corrections to the velocities if the process owns this atom
+
   if (i0 < nlocal) {
     for (int k=0; k<3; k++) 
       v[i0][k] -= imass[0] * (  l[0] * r01[k] + l[1] * r02[k] + l[2] * r03[k]);
@@ -472,13 +533,16 @@ void FixRattle::solve2x2exactly(const double a[][2],
   double determ, determinv;
 
   // calculate the determinant of the matrix
+
   determ = a[0][0] * a[1][1] - a[0][1] * a[1][0];
 
   // check if matrix is actually invertible
-  if (determ == 0.0) error->one(FLERR,"RATTLE determinant = 0.0");
+
+  if (determ == 0.0) error->one(FLERR,"Rattle determinant = 0.0");
   determinv = 1.0/determ;
 
-  // Calcualte the solution:  (l01, l02)^T = A^(-1) * c
+  // Calculate the solution:  (l01, l02)^T = A^(-1) * c
+
   l[0] = determinv * ( a[1][1] * c[0]  - a[0][1] * c[1]);
   l[1] = determinv * (-a[1][0] * c[0]  + a[0][0] * c[1]);
 }
@@ -492,14 +556,17 @@ void FixRattle::solve3x3exactly(const double a[][3],
   double determ, determinv;
   
   // calculate the determinant of the matrix
+
   determ = a[0][0]*a[1][1]*a[2][2] + a[0][1]*a[1][2]*a[2][0] + 
     a[0][2]*a[1][0]*a[2][1] - a[0][0]*a[1][2]*a[2][1] - 
     a[0][1]*a[1][0]*a[2][2] - a[0][2]*a[1][1]*a[2][0];
   
   // check if matrix is actually invertible
-  if (determ == 0.0) error->one(FLERR,"RATTLE determinant = 0.0");
+
+  if (determ == 0.0) error->one(FLERR,"Rattle determinant = 0.0");
 
   // calculate the inverse 3x3 matrix: A^(-1) = (ai_jk)
+
   determinv = 1.0/determ;
   ai[0][0] =  determinv * (a[1][1]*a[2][2] - a[1][2]*a[2][1]);
   ai[0][1] = -determinv * (a[0][1]*a[2][2] - a[0][2]*a[2][1]); 
@@ -512,6 +579,7 @@ void FixRattle::solve3x3exactly(const double a[][3],
   ai[2][2] =  determinv * (a[0][0]*a[1][1] - a[0][1]*a[1][0]);
 
   // calculate the solution:  (l01, l02, l12)^T = A^(-1) * c
+
   for (int i=0; i<3; i++) {
     l[i] = 0;
     for (int j=0; j<3; j++) 
@@ -534,6 +602,8 @@ void FixRattle::reset_dt()
 void FixRattle::update_v_half_nocons()
 {
   double dtfvinvm;
+  dtfv = 0.5 * update->dt * force->ftm2v;
+
   if (rmass) {
     for (int i = 0; i < nlocal; i++) {
       if (shake_flag[i]) {
@@ -563,9 +633,11 @@ void FixRattle::update_v_half_nocons()
 void FixRattle::update_v_half_nocons_respa(int ilevel)
 {
   // select timestep for current level
+
   dtfv = 0.5 * step_respa[ilevel] * force->ftm2v;
 
   // carry out unconstrained velocity update
+
   update_v_half_nocons();
 }
 
@@ -621,15 +693,6 @@ int FixRattle::pack_forward_comm(int n, int *list, double *buf,
         buf[m++] = v[j][2];
       }
       break;
-
-    case X:
-      for (i = 0; i < n; i++) {
-        j = list[i];
-        buf[m++] = x[j][0];
-        buf[m++] = x[j][1];
-        buf[m++] = x[j][2];
-      }
-      break;
   }
   return m;
 } 
@@ -662,25 +725,74 @@ void FixRattle::unpack_forward_comm(int n, int first, double *buf)
         v[i][2] = buf[m++];
       }
       break;
-
-    case X:
-      for (i = first; i < last; i++) {
-        x[i][0] = buf[m++];
-        x[i][1] = buf[m++];
-        x[i][2] = buf[m++];
-      }
-      break;
   }
 }
 
 
 /* ----------------------------------------------------------------------
-   wrapper method for end_of_step fixes which modify the coordinates
+  Let shake calculate new constraining forces for the coordinates;
+  As opposed to the regular shake call, this method is usually called from
+  end_of_step fixes after the second velocity integration has happened.
 ------------------------------------------------------------------------- */
 
-void FixRattle::coordinate_constraints_end_of_step() {
+void FixRattle::shake_end_of_step(int vflag) {
+
+  if (nprocs > 1) {
+    comm_mode = V;
+    comm->forward_comm_fix(this);
+  }
+
   comm_mode = XSHAKE;
-  FixShake::coordinate_constraints_end_of_step();
+  FixShake::shake_end_of_step(vflag);
+}
+
+
+/* ----------------------------------------------------------------------
+  Let shake calculate new constraining forces and correct the 
+  coordinates. Nothing to do for rattle here.
+------------------------------------------------------------------------- */
+
+void FixRattle::correct_coordinates(int vflag) {
+  comm_mode = XSHAKE;
+  FixShake::correct_coordinates(vflag);
+}
+
+
+/* ----------------------------------------------------------------------
+   Remove the velocity component along any bond.
+------------------------------------------------------------------------- */
+
+void FixRattle::correct_velocities() {
+
+  // Copy current velocities instead of unconstrained_update, because the correction
+  // should happen instantaneously and not after the next half step.
+
+  for (int i = 0; i < atom->nlocal; i++) {
+    if (shake_flag[i]) {
+      for (int k=0; k<3; k++)
+        vp[i][k] = v[i][k];
+    }
+    else
+      vp[i][0] = vp[i][1] = vp[i][2] = 0;
+  }
+
+  // communicate the unconstrained velocities
+
+  if (nprocs > 1) {
+    comm_mode = VP;
+    comm->forward_comm_fix(this);
+  }
+
+  // correct the velocity for each molecule accordingly
+
+  int m;
+  for (int i = 0; i < nlist; i++) {
+    m = list[i];
+    if      (shake_flag[m] == 2)        vrattle2(m);
+    else if (shake_flag[m] == 3)        vrattle3(m);
+    else if (shake_flag[m] == 4)        vrattle4(m);
+    else                                vrattle3angle(m);
+  }
 }
 
 
@@ -695,18 +807,16 @@ void FixRattle::coordinate_constraints_end_of_step() {
 
 void FixRattle::end_of_step()
 {
-  // communicate velocities and coordinates (x and v)
   if (nprocs > 1) {
        comm_mode = V;
        comm->forward_comm_fix(this);
-    comm_mode = X;
-    comm->forward_comm_fix(this);
   }
-  if (!check_constraints(v, RATTLE_TEST_POS, RATTLE_TEST_VEL)) {
+  if (!check_constraints(v, RATTLE_TEST_POS, RATTLE_TEST_VEL) && RATTLE_RAISE_ERROR) {
-    error->one(FLERR, "RATTLE failed");
+    error->one(FLERR, "Rattle failed ");
   }
 }
 
+
 /* ---------------------------------------------------------------------- */
 
 bool FixRattle::check_constraints(double **v, bool checkr, bool checkv)
@@ -721,6 +831,7 @@ bool FixRattle::check_constraints(double **v, bool checkr, bool checkv)
     else if (shake_flag[m] == 4)     ret =   check4(v,m,checkr,checkv); 
     else                             ret =   check3angle(v,m,checkr,checkv);
     i++;
+    if (!RATTLE_RAISE_ERROR)         ret = true;
   }
   return ret;
 }
@@ -743,13 +854,13 @@ bool FixRattle::check2(double **v, int m, bool checkr, bool checkv)
 
   stat = !(checkr && (fabs(sqrt(MathExtra::dot3(r01,r01)) - bond1) > tol));
   if (!stat)   
-     error->one(FLERR,"RATTLE coordinate constraints are not satisfied "
+     error->one(FLERR,"Coordinate constraints are not satisfied "
-                "up to desired tolerance");
+                "up to desired tolerance ");
 
   stat = !(checkv && (fabs(MathExtra::dot3(r01,v01)) > tol));
   if (!stat)   
-     error->one(FLERR,"RATTLE velocity constraints are not satisfied "
+     error->one(FLERR,"Velocity constraints are not satisfied "
-                "up to desired tolerance");
+                "up to desired tolerance ");
   return stat;
 }
 
@@ -780,14 +891,14 @@ bool FixRattle::check3(double **v, int m, bool checkr, bool checkv)
   stat = !(checkr && (fabs(sqrt(MathExtra::dot3(r01,r01)) - bond1) > tol ||
                       fabs(sqrt(MathExtra::dot3(r02,r02))-bond2) > tol));
   if (!stat)   
-     error->one(FLERR,"RATTLE coordinate constraints are not satisfied "
+     error->one(FLERR,"Coordinate constraints are not satisfied "
-                "up to desired tolerance");
+                "up to desired tolerance ");
 
   stat = !(checkv && (fabs(MathExtra::dot3(r01,v01)) > tol || 
                       fabs(MathExtra::dot3(r02,v02)) > tol));
   if (!stat)   
-     error->one(FLERR,"RATTLE velocity constraints are not satisfied "
+     error->one(FLERR,"Velocity constraints are not satisfied "
-                "up to desired tolerance");
+                "up to desired tolerance ");
   return stat;
 }
 
@@ -823,15 +934,15 @@ bool FixRattle::check4(double **v, int m, bool checkr, bool checkv)
                       fabs(sqrt(MathExtra::dot3(r02,r02))-bond2) > tol || 
                       fabs(sqrt(MathExtra::dot3(r03,r03))-bond3) > tol));
   if (!stat)   
-     error->one(FLERR,"RATTLE coordinate constraints are not satisfied "
+     error->one(FLERR,"Coordinate constraints are not satisfied "
-                "up to desired tolerance");
+                "up to desired tolerance ");
 
   stat = !(checkv && (fabs(MathExtra::dot3(r01,v01)) > tol || 
                       fabs(MathExtra::dot3(r02,v02)) > tol || 
                       fabs(MathExtra::dot3(r03,v03)) > tol));
   if (!stat)   
-     error->one(FLERR,"RATTLE velocity constraints are not satisfied "
+     error->one(FLERR,"Velocity constraints are not satisfied "
-                "up to desired tolerance");
+                "up to desired tolerance ");
   return stat;
 }
 
@@ -862,18 +973,44 @@ bool FixRattle::check3angle(double **v, int m, bool checkr, bool checkv)
   MathExtra::sub3(v[i2],v[i0],v02);
   MathExtra::sub3(v[i2],v[i1],v12);
 
-  stat = !(checkr && (fabs(sqrt(MathExtra::dot3(r01,r01)) - bond1) > tol ||
-                      fabs(sqrt(MathExtra::dot3(r02,r02))-bond2) > tol ||
-                      fabs(sqrt(MathExtra::dot3(r12,r12))-bond12) > tol));
-  if (!stat)
-     error->one(FLERR,"RATTLE coordinate constraints are not satisfied "
-                "up to desired tolerance");
 
-  stat = !(checkv && (fabs(MathExtra::dot3(r01,v01)) > tol ||
+
-                      fabs(MathExtra::dot3(r02,v02)) > tol ||
+  double db1 = fabs(sqrt(MathExtra::dot3(r01,r01)) - bond1);
-                      fabs(MathExtra::dot3(r12,v12)) > tol));
+  double db2 = fabs(sqrt(MathExtra::dot3(r02,r02))-bond2);
-  if (!stat)
+  double db12 = fabs(sqrt(MathExtra::dot3(r12,r12))-bond12);
-     error->one(FLERR,"RATTLE velocity constraints are not satisfied "
+
-                "up to desired tolerance");
+
+  stat = !(checkr && (db1 > tol || 
+                      db2 > tol || 
+                      db12 > tol));
+
+  if (derr_max < db1/bond1)    derr_max = db1/bond1;
+  if (derr_max < db2/bond2)    derr_max = db2/bond2;
+  if (derr_max < db12/bond12)  derr_max = db12/bond12;
+
+
+  if (!stat && RATTLE_RAISE_ERROR)   
+     error->one(FLERR,"Coordinate constraints are not satisfied "
+                "up to desired tolerance ");
+
+  double dv1 = fabs(MathExtra::dot3(r01,v01));
+  double dv2 = fabs(MathExtra::dot3(r02,v02));
+  double dv12 = fabs(MathExtra::dot3(r12,v12));
+
+  if (verr_max < dv1)    verr_max = dv1;
+  if (verr_max < dv2)    verr_max = dv2;
+  if (verr_max < dv12)   verr_max = dv12;
+
+
+  stat = !(checkv && (dv1 > tol || 
+                      dv2 > tol || 
+                      dv12> tol));
+
+
+  if (!stat && RATTLE_RAISE_ERROR)   
+     error->one(FLERR,"Velocity constraints are not satisfied "
+                "up to desired tolerance!");
+
+
   return stat;
 }

src/RIGID/fix_rattle.h

@@ -30,22 +30,27 @@ class FixRattle : public FixShake {
   double **vp;                // array for unconstrained velocities
   double dtfv;                // timestep for velocity update
   int comm_mode;              // mode for communication pack/unpack
+  double derr_max;            // distance error
+  double verr_max;            // velocity error
 
   FixRattle(class LAMMPS *, int, char **);
   ~FixRattle();
   int setmask();
-  void init();
+  virtual void init();
-  void post_force(int);
+  virtual void post_force(int);
-  void post_force_respa(int, int, int);
+  virtual void post_force_respa(int, int, int);
-  void final_integrate();
+  virtual void final_integrate();
-  void final_integrate_respa(int,int);
+  virtual void final_integrate_respa(int,int);
-  void coordinate_constraints_end_of_step();
 
-  double memory_usage();
+  virtual void correct_coordinates(int vflag);
-  void grow_arrays(int);
+  virtual void correct_velocities();
-  int pack_forward_comm(int, int *, double *, int, int *);
+  virtual void shake_end_of_step(int vflag);
-  void unpack_forward_comm(int, int, double *);
+
-  void reset_dt();
+  virtual double memory_usage();
+  virtual void grow_arrays(int);
+  virtual int pack_forward_comm(int, int *, double *, int, int *);
+  virtual void unpack_forward_comm(int, int, double *);
+  virtual void reset_dt();
 
  private:
   void update_v_half_nocons();
@@ -58,7 +63,7 @@ class FixRattle : public FixShake {
   void solve3x3exactly(const double a[][3], const double c[], double l[]);
   void solve2x2exactly(const double a[][2], const double c[], double l[]);
 
-  // debugging methosd
+  // debugging methods
 
   bool check3angle(double ** v, int m, bool checkr, bool checkv);
   bool check2(double **v, int m, bool checkr, bool checkv);
@@ -73,6 +78,7 @@ class FixRattle : public FixShake {
 #endif
 #endif
 
+
 /* ERROR/WARNING messages:
 
 W: Fix rattle should come after all other integration fixes
@@ -82,20 +88,20 @@ atom positions.  Thus it should be the last integration fix specified.
 If not, it will not satisfy the desired constraints as well as it
 otherwise would.
 
-E: RATTLE determinant = 0.0
+E: Rattle determinant = 0.0
 
 The determinant of the matrix being solved for a single cluster
 specified by the fix rattle command is numerically invalid.
 
-E: RATTLE failed
+E: Rattle failed
 
 Certain constraints were not satisfied.
 
-E: RATTLE coordinate constraints are not satisfied up to desired tolerance
+E: Coordinate constraints are not satisfied up to desired tolerance
 
 Self-explanatory.
 
-E: RATTLE velocity constraints are not satisfied up to desired tolerance
+E: Rattle velocity constraints are not satisfied up to desired tolerance
 
 Self-explanatory.
 

src/RIGID/fix_shake.cpp

@@ -17,6 +17,7 @@
 #include <string.h>
 #include <stdio.h>
 #include "fix_shake.h"
+#include "fix_rattle.h"
 #include "atom.h"
 #include "atom_vec.h"
 #include "molecule.h"
@@ -56,7 +57,6 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
   virial_flag = 1;
   create_attribute = 1;
   dof_flag = 1;
-
   // error check
 
   molecular = atom->molecular;
@@ -71,6 +71,9 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
   shake_type = NULL;
   xshake = NULL;
 
+  ftmp = NULL; 
+  vtmp = NULL;
+
   grow_arrays(atom->nmax);
   atom->add_callback(0);
 
@@ -206,7 +209,6 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
   // SHAKE vs RATTLE
 
   rattle = 0;
-  vflag_post_force = 0;
 
   // identify all SHAKE clusters
 
@@ -255,6 +257,9 @@ FixShake::~FixShake()
   memory->destroy(shake_atom);
   memory->destroy(shake_type);
   memory->destroy(xshake);
+  memory->destroy(ftmp);
+  memory->destroy(vtmp);
+
 
   delete [] bond_flag;
   delete [] angle_flag;
@@ -433,30 +438,27 @@ void FixShake::setup(int vflag)
       next_output = (ntimestep/output_every)*output_every + output_every;
   } else next_output = -1;
 
-  // half timestep constraint on pre-step, full timestep thereafter
 
-  if (strstr(update->integrate_style,"verlet")) {
+  // set respa to 0 if verlet is used and to 1 otherwise
+
+  if (strstr(update->integrate_style,"verlet")) 
     respa = 0;
-    dtv     = update->dt;
+  else
-    dtfsq   = 0.5 * update->dt * update->dt * force->ftm2v;
-    FixShake::post_force(vflag);
-    if (!rattle) dtfsq = update->dt * update->dt * force->ftm2v;
-
-  } else {
     respa = 1;
-    dtv = step_respa[0];
-    dtf_innerhalf = 0.5 * step_respa[0] * force->ftm2v;
-    dtf_inner = dtf_innerhalf;
 
-    // apply correction to all rRESPA levels
 
-    for (int ilevel = 0; ilevel < nlevels_respa; ilevel++) {
+  // correct geometry of cluster if necessary
-      ((Respa *) update->integrate)->copy_flevel_f(ilevel);
+
-      FixShake::post_force_respa(vflag,ilevel,loop_respa[ilevel]-1);
+  correct_coordinates(vflag);
-      ((Respa *) update->integrate)->copy_f_flevel(ilevel);
+
-    }
+  // remove velocities along any bonds
-    if (!rattle) dtf_inner = step_respa[0] * force->ftm2v;
+
-  }
+  correct_velocities();
+
+  // precalculate constraining forces for first integration step
+
+  shake_end_of_step(vflag);
+
 }
 
 /* ----------------------------------------------------------------------
@@ -572,7 +574,6 @@ void FixShake::post_force(int vflag)
   }
   
   // store vflag for coordinate_constraints_end_of_step()
-
   vflag_post_force = vflag;
 }
 
@@ -620,7 +621,6 @@ void FixShake::post_force_respa(int vflag, int ilevel, int iloop)
   }
 
   // store vflag for coordinate_constraints_end_of_step()
-
   vflag_post_force = vflag;
 }
 
@@ -2397,6 +2397,11 @@ void FixShake::grow_arrays(int nmax)
   memory->grow(shake_type,nmax,3,"shake:shake_type");
   memory->destroy(xshake);
   memory->create(xshake,nmax,3,"shake:xshake");
+  memory->destroy(ftmp);
+  memory->create(ftmp,nmax,3,"shake:ftmp");
+  memory->destroy(vtmp);
+  memory->create(vtmp,nmax,3,"shake:vtmp");
+
 }
 
 /* ----------------------------------------------------------------------
@@ -2683,25 +2688,116 @@ void *FixShake::extract(const char *str, int &dim)
 
 
 /* ----------------------------------------------------------------------
-   wrapper method for end_of_step fixes which modify the coordinates
+   Add coordinate constraining forces; this method is called
+   at the end of a timestep.
 ------------------------------------------------------------------------- */
 
-void FixShake::coordinate_constraints_end_of_step() {
+void FixShake::shake_end_of_step(int vflag) {
+
   if (!respa) {
+    dtv     = update->dt;
     dtfsq   = 0.5 * update->dt * update->dt * force->ftm2v;
-    FixShake::post_force(vflag_post_force);
+    FixShake::post_force(vflag);
     if (!rattle) dtfsq = update->dt * update->dt * force->ftm2v;
-  }
+
-  else {
+  } else {
+    dtv = step_respa[0];
     dtf_innerhalf = 0.5 * step_respa[0] * force->ftm2v;
     dtf_inner = dtf_innerhalf;
+
     // apply correction to all rRESPA levels
+
     for (int ilevel = 0; ilevel < nlevels_respa; ilevel++) {
       ((Respa *) update->integrate)->copy_flevel_f(ilevel);
-      FixShake::post_force_respa(vflag_post_force,ilevel,loop_respa[ilevel]-1);
+      FixShake::post_force_respa(vflag,ilevel,loop_respa[ilevel]-1);
       ((Respa *) update->integrate)->copy_f_flevel(ilevel);
     }
     if (!rattle) dtf_inner = step_respa[0] * force->ftm2v;
   }
 }
 
+
+
+/* ----------------------------------------------------------------------
+   wrapper method for end_of_step fixes which modify velocities
+------------------------------------------------------------------------- */
+
+void FixShake::correct_velocities() { 
+
+}
+
+
+/* ----------------------------------------------------------------------
+   Calculate constraining forces based on the current configuration
+   and change coordinates.
+------------------------------------------------------------------------- */
+
+void FixShake::correct_coordinates(int vflag) { 
+   
+  // save current forces and velocities so that you 
+  // initialise them to zero such that FixShake::unconstrained_coordinate_update has no effect 
+
+  for (int j=0; j<nlocal; j++) {
+    for (int k=0; k<3; k++) {
+
+      // store current value of forces and velocities
+
+      ftmp[j][k] = f[j][k];
+      vtmp[j][k] = v[j][k];
+
+      // set f and v to zero for SHAKE
+
+      v[j][k] = 0;
+      f[j][k] = 0;
+    }
+  }
+
+  // call SHAKE to correct the coordinates which were updated without constraints
+  // IMPORTANT: use 1 as argument and thereby enforce velocity Verlet
+
+  dtfsq   = 0.5 * update->dt * update->dt * force->ftm2v;
+  FixShake::post_force(vflag);
+
+  // integrate coordiantes: x' = xnp1 + dt^2/2m_i * f, where f is the constraining force
+  // NOTE: After this command, the coordinates geometry of the molecules will be correct! 
+
+  double dtfmsq; 
+  if (rmass) {  
+    for (int i = 0; i < nlocal; i++) {  
+      dtfmsq = dtfsq/ rmass[i];  
+      x[i][0] = x[i][0] + dtfmsq*f[i][0];  
+      x[i][1] = x[i][1] + dtfmsq*f[i][1];  
+      x[i][2] = x[i][2] + dtfmsq*f[i][2];  
+    }  
+  }  
+  else {  
+    for (int i = 0; i < nlocal; i++) {  
+      dtfmsq = dtfsq / mass[type[i]];  
+      x[i][0] = x[i][0] + dtfmsq*f[i][0];  
+      x[i][1] = x[i][1] + dtfmsq*f[i][1];  
+      x[i][2] = x[i][2] + dtfmsq*f[i][2];  
+    }  
+  }
+
+  // copy forces and velocities back
+
+  for (int j=0; j<nlocal; j++) {
+    for (int k=0; k<3; k++) {
+      f[j][k] = ftmp[j][k];
+      v[j][k] = vtmp[j][k];
+    }
+  }
+
+  if (!rattle) dtfsq = update->dt * update->dt * force->ftm2v;
+
+  // communicate changes
+  // NOTE: for compatibility xshake is temporarily set to x, such that pack/unpack_forward 
+  //       can be used for communicating the coordinates. 
+
+  double **xtmp = xshake;
+  xshake = x;
+  if (nprocs > 1) {
+    comm->forward_comm_fix(this);
+  }
+  xshake = xtmp;
+}

src/RIGID/fix_shake.h

@@ -25,6 +25,9 @@ FixStyle(shake,FixShake)
 namespace LAMMPS_NS {
 
 class FixShake : public Fix {
+
+ friend class FixEHEX;
+
  public:
   FixShake(class LAMMPS *, int, char **);
   virtual ~FixShake();
@@ -46,7 +49,13 @@ class FixShake : public Fix {
   virtual int unpack_exchange(int, double *);
   virtual int pack_forward_comm(int, int *, double *, int, int *);
   virtual void unpack_forward_comm(int, int, double *);
-  virtual void coordinate_constraints_end_of_step();
+
+
+
+  virtual void shake_end_of_step(int vflag);
+  virtual void correct_coordinates(int vflag);
+  virtual void correct_velocities();
+
 
   int dof(int);
   virtual void reset_dt();
@@ -77,6 +86,8 @@ class FixShake : public Fix {
   double *step_respa;
 
   double **x,**v,**f;                    // local ptrs to atom class quantities
+  double **ftmp, **vtmp;                 // pointers to temporary arrays for forces and velocities
+
   double *mass,*rmass;
   int *type;
   int nlocal;