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

src/RIGID/fix_shake.cpp

@@ -205,6 +205,7 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
   // SHAKE vs RATTLE
 
   rattle = 0;
+  vflag_post_force = 0;
 
   // identify all SHAKE clusters
 
@@ -432,12 +433,14 @@ void FixShake::setup(int vflag)
   // half timestep constraint on pre-step, full timestep thereafter
 
   if (strstr(update->integrate_style,"verlet")) {
+    respa   = 0;   
     dtv     = update->dt;
     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;
@@ -564,6 +567,10 @@ void FixShake::post_force(int vflag)
     else if (shake_flag[m] == 4) shake4(m);
     else shake3angle(m);
   }
+  
+  // store vflag for coordinate_constraints_end_of_step()
+
+  vflag_post_force = vflag;
 }
 
 /* ----------------------------------------------------------------------
@@ -608,6 +615,10 @@ void FixShake::post_force_respa(int vflag, int ilevel, int iloop)
     else if (shake_flag[m] == 4) shake4(m);
     else shake3angle(m);
   }
+
+  // store vflag for coordinate_constraints_end_of_step()
+
+  vflag_post_force = vflag;
 }
 
 /* ----------------------------------------------------------------------
@@ -2665,3 +2676,29 @@ void *FixShake::extract(const char *str, int &dim)
   }
   return NULL;
 }
+
+
+
+/* ----------------------------------------------------------------------
+   wrapper method for end_of_step fixes which modify the coordinates
+------------------------------------------------------------------------- */
+
+void FixShake::coordinate_constraints_end_of_step() {
+  if (!respa) {
+    dtfsq   = 0.5 * update->dt * update->dt * force->ftm2v;
+    FixShake::post_force(vflag_post_force);
+    if (!rattle) dtfsq = update->dt * update->dt * force->ftm2v;
+  } 
+  else {
+    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);
+      ((Respa *) update->integrate)->copy_f_flevel(ilevel);
+    }
+    if (!rattle) dtf_inner = step_respa[0] * force->ftm2v;
+  }
+}
+

src/RIGID/fix_shake.h

@@ -46,19 +46,22 @@ 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();
 
   int dof(int);
   virtual void reset_dt();
   void *extract(const char *, int &);
 
  protected:
+  int vflag_post_force;                  // store the vflag of last post_force call
+  int respa;                             // 0 = vel. Verlet, 1 = respa
   int me,nprocs;
   int rattle;                            // 0 = SHAKE, 1 = RATTLE
   double tolerance;                      // SHAKE tolerance
   int max_iter;                          // max # of SHAKE iterations
   int output_every;                      // SHAKE stat output every so often
   bigint next_output;                    // timestep for next output
-
+  
                                          // settings from input command
   int *bond_flag,*angle_flag;            // bond/angle types to constrain
   int *type_flag;                        // constrain bonds to these types

src/fix_respa.h

@@ -50,7 +50,6 @@ class FixRespa : public Fix {
 
 #endif
 #endif
-
 /* ERROR/WARNING messages:
 
 */