diff --git a/src/RIGID/fix_shake.cpp b/src/RIGID/fix_shake.cpp index 0dd6c85d5f..bf41b17a5b 100644 --- a/src/RIGID/fix_shake.cpp +++ b/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; + } +} + diff --git a/src/RIGID/fix_shake.h b/src/RIGID/fix_shake.h index 74e09d12d8..049971ca4c 100644 --- a/src/RIGID/fix_shake.h +++ b/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 diff --git a/src/fix_respa.h b/src/fix_respa.h index d65587cdc9..7c7bd53a61 100644 --- a/src/fix_respa.h +++ b/src/fix_respa.h @@ -50,7 +50,6 @@ class FixRespa : public Fix { #endif #endif - /* ERROR/WARNING messages: */