diff --git a/src/ASPHERE/fix_npt_asphere.cpp b/src/ASPHERE/fix_npt_asphere.cpp index eb4b1f8901..9f7233bc46 100755 --- a/src/ASPHERE/fix_npt_asphere.cpp +++ b/src/ASPHERE/fix_npt_asphere.cpp @@ -160,9 +160,9 @@ void FixNPTAsphere::initial_integrate(int vflag) } } - // remap simulation box and all owned atoms by 1/2 step + // remap simulation box by 1/2 step - remap(0); + remap(); // update x by full step for atoms in group @@ -192,10 +192,10 @@ void FixNPTAsphere::initial_integrate(int vflag) } } - // remap simulation box and all owned atoms by 1/2 step + // remap simulation box by 1/2 step // redo KSpace coeffs since volume has changed - remap(0); + remap(); if (kspace_flag) force->kspace->setup(); } diff --git a/src/POEMS/fix_poems.cpp b/src/POEMS/fix_poems.cpp index a8feff5342..f4edb3f1cf 100644 --- a/src/POEMS/fix_poems.cpp +++ b/src/POEMS/fix_poems.cpp @@ -775,10 +775,8 @@ void FixPOEMS::final_integrate() /* ---------------------------------------------------------------------- */ -void FixPOEMS::initial_integrate_respa(int vflag, int ilevel, int flag) +void FixPOEMS::initial_integrate_respa(int vflag, int ilevel, int iloop) { - if (flag) return; // only used by NPT,NPH - dtv = step_respa[ilevel]; dtf = 0.5 * step_respa[ilevel] * force->ftm2v; dthalf = 0.5 * step_respa[ilevel]; @@ -796,7 +794,7 @@ void FixPOEMS::post_force_respa(int vflag, int ilevel, int iloop) /* ---------------------------------------------------------------------- */ -void FixPOEMS::final_integrate_respa(int ilevel) +void FixPOEMS::final_integrate_respa(int ilevel, int iloop) { dtf = 0.5 * step_respa[ilevel] * force->ftm2v; final_integrate(); diff --git a/src/POEMS/fix_poems.h b/src/POEMS/fix_poems.h index 2374a1626f..8d41964099 100644 --- a/src/POEMS/fix_poems.h +++ b/src/POEMS/fix_poems.h @@ -36,7 +36,7 @@ class FixPOEMS : public Fix { void final_integrate(); void initial_integrate_respa(int, int, int); void post_force_respa(int, int, int); - void final_integrate_respa(int); + void final_integrate_respa(int, int); void grow_arrays(int); void copy_arrays(int, int); diff --git a/src/atom.cpp b/src/atom.cpp index 22019674df..ecebf7a52b 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -1496,7 +1496,8 @@ void Atom::add_callback(int flag) // find the fix // if find NULL ptr: - // it's this one, since it is deleted at this point in creation + // it's this one, since it is being replaced and has just been deleted + // at this point in re-creation // if don't find NULL ptr: // i is set to nfix = new one currently being added at end of list diff --git a/src/fix.h b/src/fix.h index dcbcb08323..0698ffcb3c 100644 --- a/src/fix.h +++ b/src/fix.h @@ -115,7 +115,7 @@ class Fix : protected Pointers { virtual void post_integrate_respa(int, int) {} virtual void pre_force_respa(int, int, int) {} virtual void post_force_respa(int, int, int) {} - virtual void final_integrate_respa(int) {} + virtual void final_integrate_respa(int, int) {} virtual void min_pre_exchange() {} virtual void min_post_force(int) {} diff --git a/src/fix_move.cpp b/src/fix_move.cpp index 38a4bd4be4..3a68079ff9 100644 --- a/src/fix_move.cpp +++ b/src/fix_move.cpp @@ -769,10 +769,8 @@ void FixMove::final_integrate() /* ---------------------------------------------------------------------- */ -void FixMove::initial_integrate_respa(int vflag, int ilevel, int flag) +void FixMove::initial_integrate_respa(int vflag, int ilevel, int iloop) { - if (flag) return; // only used by NPT,NPH - // outermost level - update v and x // all other levels - nothing @@ -781,7 +779,7 @@ void FixMove::initial_integrate_respa(int vflag, int ilevel, int flag) /* ---------------------------------------------------------------------- */ -void FixMove::final_integrate_respa(int ilevel) +void FixMove::final_integrate_respa(int ilevel, int iloop) { if (ilevel == nlevels_respa-1) final_integrate(); } diff --git a/src/fix_move.h b/src/fix_move.h index 9d87916c19..2bf0370c94 100644 --- a/src/fix_move.h +++ b/src/fix_move.h @@ -33,7 +33,7 @@ class FixMove : public Fix { void initial_integrate(int); void final_integrate(); void initial_integrate_respa(int, int, int); - void final_integrate_respa(int); + void final_integrate_respa(int, int); double memory_usage(); void grow_arrays(int); diff --git a/src/fix_nph.cpp b/src/fix_nph.cpp index 39b5c0c6ea..5169ce2943 100644 --- a/src/fix_nph.cpp +++ b/src/fix_nph.cpp @@ -407,9 +407,9 @@ void FixNPH::initial_integrate(int vflag) } } - // remap simulation box and all owned atoms by 1/2 step + // remap simulation box by 1/2 step - remap(0); + remap(); // x update by full step only for atoms in NPH group @@ -421,10 +421,10 @@ void FixNPH::initial_integrate(int vflag) } } - // remap simulation box and all owned atoms by 1/2 step + // remap simulation box by 1/2 step // redo KSpace coeffs since volume has changed - remap(0); + remap(); if (kspace_flag) force->kspace->setup(); } @@ -499,18 +499,8 @@ void FixNPH::final_integrate() /* ---------------------------------------------------------------------- */ -void FixNPH::initial_integrate_respa(int vflag, int ilevel, int flag) +void FixNPH::initial_integrate_respa(int vflag, int ilevel, int iloop) { - // if flag = 1, then is 2nd call at outermost level from rRESPA - // perform 2nd half of box remap on own + ghost atoms and return - // redo KSpace coeffs since volume has changed - - if (flag == 1) { - remap(1); - if (kspace_flag) force->kspace->setup(); - return; - } - // set timesteps by level double dtfm; @@ -578,9 +568,10 @@ void FixNPH::initial_integrate_respa(int vflag, int ilevel, int flag) } } - // remap simulation box and all owned atoms by 1/2 step + // remap simulation box by 1/2 step - remap(0); + remap(); + remap2flag = 1; } else { @@ -607,7 +598,11 @@ void FixNPH::initial_integrate_respa(int vflag, int ilevel, int flag) } } - // innermost level - also update x only for atoms in NPH group + // innermost level - also update x for atoms in group + // if remap2flag: + // this is 1st call at innermost level from rRESPA after 1st half remap + // perform 2nd half of box remap + // redo KSpace coeffs since volume has changed if (ilevel == 0) { for (int i = 0; i < nlocal; i++) { @@ -617,12 +612,17 @@ void FixNPH::initial_integrate_respa(int vflag, int ilevel, int flag) x[i][2] += dtv * v[i][2]; } } + if (remap2flag) { + remap(); + if (kspace_flag) force->kspace->setup(); + remap2flag = 0; + } } } /* ---------------------------------------------------------------------- */ -void FixNPH::final_integrate_respa(int ilevel) +void FixNPH::final_integrate_respa(int ilevel, int iloop) { double dtfm; @@ -699,26 +699,24 @@ void FixNPH::couple() /* ---------------------------------------------------------------------- change box size - remap owned or owned+ghost atoms depending on flag remap all atoms or fix group atoms depending on allremap flag if rigid bodies exist, scale rigid body centers-of-mass ------------------------------------------------------------------------- */ -void FixNPH::remap(int flag) +void FixNPH::remap() { - int i,n; + int i; double oldlo,oldhi,ctr; double **x = atom->x; int *mask = atom->mask; - if (flag) n = atom->nlocal + atom->nghost; - else n = atom->nlocal; + int nlocal = atom->nlocal; // convert pertinent atoms and rigid bodies to lamda coords - if (allremap) domain->x2lamda(n); + if (allremap) domain->x2lamda(nlocal); else { - for (i = 0; i < n; i++) + for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) domain->x2lamda(x[i],x[i]); } @@ -744,9 +742,9 @@ void FixNPH::remap(int flag) // convert pertinent atoms and rigid bodies back to box coords - if (allremap) domain->lamda2x(n); + if (allremap) domain->lamda2x(nlocal); else { - for (i = 0; i < n; i++) + for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) domain->lamda2x(x[i],x[i]); } diff --git a/src/fix_nph.h b/src/fix_nph.h index 78dbe04390..1911e1d739 100644 --- a/src/fix_nph.h +++ b/src/fix_nph.h @@ -34,7 +34,7 @@ class FixNPH : public Fix { void initial_integrate(int); void final_integrate(); void initial_integrate_respa(int, int, int); - void final_integrate_respa(int); + void final_integrate_respa(int, int); double compute_scalar(); void write_restart(FILE *); void restart(char *); @@ -60,13 +60,14 @@ class FixNPH : public Fix { int nlevels_respa; double *step_respa; + int remap2flag; // flag for performing 2nd half remap() char *id_temp,*id_press; class Compute *temperature,*pressure; int tflag,pflag; void couple(); - void remap(int); + void remap(); }; } diff --git a/src/fix_npt.cpp b/src/fix_npt.cpp index a188aec11d..bbd534ee8e 100644 --- a/src/fix_npt.cpp +++ b/src/fix_npt.cpp @@ -453,9 +453,9 @@ void FixNPT::initial_integrate(int vflag) } } - // remap simulation box and all owned atoms by 1/2 step + // remap simulation box by 1/2 step - remap(0); + remap(); // update x by full step for atoms in group @@ -467,10 +467,10 @@ void FixNPT::initial_integrate(int vflag) } } - // remap simulation box and all owned atoms by 1/2 step + // remap simulation box by 1/2 step // redo KSpace coeffs since volume has changed - remap(0); + remap(); if (kspace_flag) force->kspace->setup(); } @@ -585,18 +585,8 @@ void FixNPT::final_integrate() /* ---------------------------------------------------------------------- */ -void FixNPT::initial_integrate_respa(int vflag, int ilevel, int flag) +void FixNPT::initial_integrate_respa(int vflag, int ilevel, int iloop) { - // if flag = 1, then is 2nd call at outermost level from rRESPA - // perform 2nd half of box remap on own + ghost atoms and return - // redo KSpace coeffs since volume has changed - - if (flag == 1) { - remap(1); - if (kspace_flag) force->kspace->setup(); - return; - } - // set timesteps by level double dtfm; @@ -673,9 +663,10 @@ void FixNPT::initial_integrate_respa(int vflag, int ilevel, int flag) } } - // remap simulation box and all owned atoms by 1/2 step + // remap simulation box by 1/2 step - remap(0); + remap(); + remap2flag = 1; } else { @@ -704,6 +695,10 @@ void FixNPT::initial_integrate_respa(int vflag, int ilevel, int flag) } // innermost level - also update x for atoms in group + // if remap2flag: + // this is 1st call at innermost level from rRESPA after 1st half remap + // perform 2nd half of box remap + // redo KSpace coeffs since volume has changed if (ilevel == 0) { for (int i = 0; i < nlocal; i++) { @@ -713,12 +708,17 @@ void FixNPT::initial_integrate_respa(int vflag, int ilevel, int flag) x[i][2] += dtv * v[i][2]; } } + if (remap2flag) { + remap(); + if (kspace_flag) force->kspace->setup(); + remap2flag = 0; + } } } /* ---------------------------------------------------------------------- */ -void FixNPT::final_integrate_respa(int ilevel) +void FixNPT::final_integrate_respa(int ilevel, int iloop) { double dtfm; @@ -798,26 +798,24 @@ void FixNPT::couple() /* ---------------------------------------------------------------------- change box size - remap owned or owned+ghost atoms depending on flag remap all atoms or fix group atoms depending on allremap flag if rigid bodies exist, scale rigid body centers-of-mass ------------------------------------------------------------------------- */ -void FixNPT::remap(int flag) +void FixNPT::remap() { - int i,n; + int i; double oldlo,oldhi,ctr; double **x = atom->x; int *mask = atom->mask; - if (flag) n = atom->nlocal + atom->nghost; - else n = atom->nlocal; + int nlocal = atom->nlocal; // convert pertinent atoms and rigid bodies to lamda coords - if (allremap) domain->x2lamda(n); + if (allremap) domain->x2lamda(nlocal); else { - for (i = 0; i < n; i++) + for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) domain->x2lamda(x[i],x[i]); } @@ -843,9 +841,9 @@ void FixNPT::remap(int flag) // convert pertinent atoms and rigid bodies back to box coords - if (allremap) domain->lamda2x(n); + if (allremap) domain->lamda2x(nlocal); else { - for (i = 0; i < n; i++) + for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) domain->lamda2x(x[i],x[i]); } diff --git a/src/fix_npt.h b/src/fix_npt.h index 6b0273fbf2..d46619d77f 100644 --- a/src/fix_npt.h +++ b/src/fix_npt.h @@ -34,7 +34,7 @@ class FixNPT : public Fix { virtual void initial_integrate(int); virtual void final_integrate(); void initial_integrate_respa(int, int, int); - void final_integrate_respa(int); + void final_integrate_respa(int, int); double compute_scalar(); void write_restart(FILE *); void restart(char *); @@ -66,13 +66,14 @@ class FixNPT : public Fix { int nlevels_respa; double *step_respa; + int remap2flag; // flag for performing 2nd half remap() char *id_temp,*id_press; class Compute *temperature,*pressure; int tflag,pflag; void couple(); - void remap(int); + void remap(); }; } diff --git a/src/fix_npt_sphere.cpp b/src/fix_npt_sphere.cpp index 3ba3bf9861..9bd08f2c40 100644 --- a/src/fix_npt_sphere.cpp +++ b/src/fix_npt_sphere.cpp @@ -191,9 +191,9 @@ void FixNPTSphere::initial_integrate(int vflag) } } - // remap simulation box and all owned atoms by 1/2 step + // remap simulation box by 1/2 step - remap(0); + remap(); // update x by full step for atoms in group @@ -259,10 +259,10 @@ void FixNPTSphere::initial_integrate(int vflag) } } - // remap simulation box and all owned atoms by 1/2 step + // remap simulation box by 1/2 step // redo KSpace coeffs since volume has changed - remap(0); + remap(); if (kspace_flag) force->kspace->setup(); } diff --git a/src/fix_nve.cpp b/src/fix_nve.cpp index cddb2f86c3..78236b2f29 100644 --- a/src/fix_nve.cpp +++ b/src/fix_nve.cpp @@ -145,10 +145,8 @@ void FixNVE::final_integrate() /* ---------------------------------------------------------------------- */ -void FixNVE::initial_integrate_respa(int vflag, int ilevel, int flag) +void FixNVE::initial_integrate_respa(int vflag, int ilevel, int iloop) { - if (flag) return; // only used by NPT,NPH - dtv = step_respa[ilevel]; dtf = 0.5 * step_respa[ilevel] * force->ftm2v; @@ -161,7 +159,7 @@ void FixNVE::initial_integrate_respa(int vflag, int ilevel, int flag) /* ---------------------------------------------------------------------- */ -void FixNVE::final_integrate_respa(int ilevel) +void FixNVE::final_integrate_respa(int ilevel, int iloop) { dtf = 0.5 * step_respa[ilevel] * force->ftm2v; final_integrate(); diff --git a/src/fix_nve.h b/src/fix_nve.h index 6f8633733f..f894b8dc6c 100644 --- a/src/fix_nve.h +++ b/src/fix_nve.h @@ -32,7 +32,7 @@ class FixNVE : public Fix { virtual void initial_integrate(int); virtual void final_integrate(); void initial_integrate_respa(int, int, int); - void final_integrate_respa(int); + void final_integrate_respa(int, int); void reset_dt(); protected: diff --git a/src/fix_nve_limit.cpp b/src/fix_nve_limit.cpp index d6ecd4860d..ad14aa1fc9 100644 --- a/src/fix_nve_limit.cpp +++ b/src/fix_nve_limit.cpp @@ -189,10 +189,8 @@ void FixNVELimit::final_integrate() /* ---------------------------------------------------------------------- */ -void FixNVELimit::initial_integrate_respa(int vflag, int ilevel, int flag) +void FixNVELimit::initial_integrate_respa(int vflag, int ilevel, int iloop) { - if (flag) return; // only used by NPT,NPH - dtv = step_respa[ilevel]; dtf = 0.5 * step_respa[ilevel] * force->ftm2v; @@ -202,7 +200,7 @@ void FixNVELimit::initial_integrate_respa(int vflag, int ilevel, int flag) /* ---------------------------------------------------------------------- */ -void FixNVELimit::final_integrate_respa(int ilevel) +void FixNVELimit::final_integrate_respa(int ilevel, int iloop) { dtf = 0.5 * step_respa[ilevel] * force->ftm2v; final_integrate(); diff --git a/src/fix_nve_limit.h b/src/fix_nve_limit.h index 3b5613fd5d..475b2e6419 100644 --- a/src/fix_nve_limit.h +++ b/src/fix_nve_limit.h @@ -32,7 +32,7 @@ class FixNVELimit : public Fix { void initial_integrate(int); void final_integrate(); void initial_integrate_respa(int, int, int); - void final_integrate_respa(int); + void final_integrate_respa(int, int); void reset_dt(); double compute_scalar(); diff --git a/src/fix_nvt.cpp b/src/fix_nvt.cpp index 74bfc1f5b3..18a3fd12c0 100644 --- a/src/fix_nvt.cpp +++ b/src/fix_nvt.cpp @@ -442,10 +442,8 @@ void FixNVT::final_integrate() /* ---------------------------------------------------------------------- */ -void FixNVT::initial_integrate_respa(int vflag, int ilevel, int flag) +void FixNVT::initial_integrate_respa(int vflag, int ilevel, int iloop) { - if (flag) return; // only used by NPT,NPH - // set timesteps by level double dtfm; @@ -567,7 +565,7 @@ void FixNVT::initial_integrate_respa(int vflag, int ilevel, int flag) /* ---------------------------------------------------------------------- */ -void FixNVT::final_integrate_respa(int ilevel) +void FixNVT::final_integrate_respa(int ilevel, int iloop) { // set timesteps by level diff --git a/src/fix_nvt.h b/src/fix_nvt.h index 95c43126c5..e7bfbfe784 100644 --- a/src/fix_nvt.h +++ b/src/fix_nvt.h @@ -34,7 +34,7 @@ class FixNVT : public Fix { virtual void initial_integrate(int); virtual void final_integrate(); virtual void initial_integrate_respa(int, int, int); - void final_integrate_respa(int); + void final_integrate_respa(int, int); double compute_scalar(); void write_restart(FILE *); void restart(char *); diff --git a/src/fix_rigid.cpp b/src/fix_rigid.cpp index d67e571ea5..051cdb46cd 100644 --- a/src/fix_rigid.cpp +++ b/src/fix_rigid.cpp @@ -1357,10 +1357,8 @@ void FixRigid::no_squish_rotate(int k, double *p, double *q, /* ---------------------------------------------------------------------- */ -void FixRigid::initial_integrate_respa(int vflag, int ilevel, int flag) +void FixRigid::initial_integrate_respa(int vflag, int ilevel, int iloop) { - if (flag) return; // only used by NPT,NPH - dtv = step_respa[ilevel]; dtf = 0.5 * step_respa[ilevel] * force->ftm2v; dtq = 0.5 * step_respa[ilevel]; @@ -1371,7 +1369,7 @@ void FixRigid::initial_integrate_respa(int vflag, int ilevel, int flag) /* ---------------------------------------------------------------------- */ -void FixRigid::final_integrate_respa(int ilevel) +void FixRigid::final_integrate_respa(int ilevel, int iloop) { dtf = 0.5 * step_respa[ilevel] * force->ftm2v; final_integrate(); diff --git a/src/fix_rigid.h b/src/fix_rigid.h index 869fb49f35..149da91c03 100644 --- a/src/fix_rigid.h +++ b/src/fix_rigid.h @@ -34,7 +34,7 @@ class FixRigid : public Fix { virtual void initial_integrate(int); virtual void final_integrate(); void initial_integrate_respa(int, int, int); - void final_integrate_respa(int); + void final_integrate_respa(int, int); double memory_usage(); void grow_arrays(int); diff --git a/src/fix_shake.cpp b/src/fix_shake.cpp index a7280f44dc..20ef80fe35 100644 --- a/src/fix_shake.cpp +++ b/src/fix_shake.cpp @@ -498,11 +498,9 @@ void FixShake::post_force(int vflag) if (update->ntimestep == next_output) stats(); // xshake = unconstrained move with current v,f - - unconstrained_update(); - // communicate results if necessary + unconstrained_update(); if (nprocs > 1) comm->forward_comm_fix(this); // virial setup @@ -510,7 +508,48 @@ void FixShake::post_force(int vflag) if (vflag) v_setup(vflag); else evflag = 0; - // loop over clusters + // loop over clusters to add constraint forces + + int m; + for (int i = 0; i < nlist; i++) { + m = list[i]; + if (shake_flag[m] == 2) shake2(m); + else if (shake_flag[m] == 3) shake3(m); + else if (shake_flag[m] == 4) shake4(m); + else shake3angle(m); + } +} + +/* ---------------------------------------------------------------------- + enforce SHAKE constraints from rRESPA + xshake prediction portion is different than Verlet +------------------------------------------------------------------------- */ + +void FixShake::post_force_respa(int vflag, int ilevel, int iloop) +{ + // call stats only on outermost level + + printf("SHAKE POST FORCE ilevel %d iloop %d\n",ilevel,iloop); + stats(); + //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 + + if (ilevel < nlevels_respa-1 && iloop == loop_respa[ilevel]-1) return; + + // xshake = unconstrained move with current v,f as function of level + // communicate results if necessary + + unconstrained_update_respa(ilevel); + if (nprocs > 1) comm->forward_comm_fix(this); + + // virial setup, only need to compute on outermost level + + if (ilevel == nlevels_respa-1 && vflag) v_setup(vflag); + else evflag = 0; + + // loop over clusters to add constraint forces int m; for (int i = 0; i < nlist; i++) { @@ -1197,6 +1236,63 @@ void FixShake::unconstrained_update() } } +/* ---------------------------------------------------------------------- + update the unconstrained position of each atom in a rRESPA step + only for SHAKE clusters, else set to 0.0 + assumes NVE update, seems to be accurate enough for NVT,NPT,NPH as well +------------------------------------------------------------------------- */ + +void FixShake::unconstrained_update_respa(int ilevel) +{ + // xshake = atom coords after next x update in innermost loop + // depends on rRESPA level + // for levels > 0 this includes more than one velocity update + // xshake = predicted position from call to this routine at level N = + // x + dt0 (v + dtN/m fN + 1/2 dt(N-1)/m f(N-1) + ... + 1/2 dt0/m f0) + // also set dtfsq = dt0*dtN so that shake2,shake3,etc can use it + + double ***f_level = ((FixRespa *) modify->fix[ifix_respa])->f_level; + dtfsq = dtf_inner * step_respa[ilevel]; + + double invmass,dtfmsq; + int jlevel; + + if (rmass) { + for (int i = 0; i < nlocal; i++) { + if (shake_flag[i]) { + invmass = 1.0 / rmass[i]; + dtfmsq = dtfsq * invmass; + xshake[i][0] = x[i][0] + dtv*v[i][0] + dtfmsq*f[i][0]; + xshake[i][1] = x[i][1] + dtv*v[i][1] + dtfmsq*f[i][1]; + xshake[i][2] = x[i][2] + dtv*v[i][2] + dtfmsq*f[i][2]; + for (jlevel = 0; jlevel < ilevel; jlevel++) { + dtfmsq = dtf_innerhalf * step_respa[jlevel] * invmass; + xshake[i][0] += dtfmsq*f_level[i][jlevel][0]; + xshake[i][1] += dtfmsq*f_level[i][jlevel][1]; + xshake[i][2] += dtfmsq*f_level[i][jlevel][2]; + } + } else xshake[i][2] = xshake[i][1] = xshake[i][0] = 0.0; + } + + } else { + for (int i = 0; i < nlocal; i++) { + if (shake_flag[i]) { + invmass = 1.0 / mass[type[i]]; + dtfmsq = dtfsq * invmass; + xshake[i][0] = x[i][0] + dtv*v[i][0] + dtfmsq*f[i][0]; + xshake[i][1] = x[i][1] + dtv*v[i][1] + dtfmsq*f[i][1]; + xshake[i][2] = x[i][2] + dtv*v[i][2] + dtfmsq*f[i][2]; + for (jlevel = 0; jlevel < ilevel; jlevel++) { + dtfmsq = dtf_innerhalf * step_respa[jlevel] * invmass; + xshake[i][0] += dtfmsq*f_level[i][jlevel][0]; + xshake[i][1] += dtfmsq*f_level[i][jlevel][1]; + xshake[i][2] += dtfmsq*f_level[i][jlevel][2]; + } + } else xshake[i][2] = xshake[i][1] = xshake[i][0] = 0.0; + } + } +} + /* ---------------------------------------------------------------------- */ void FixShake::shake2(int m) @@ -1935,6 +2031,10 @@ void FixShake::shake3angle(int m) v_tally(nlist,list,3.0,v); } + + if (i0 < 20) + printf("AAA %d %d %d %g %g %g: %g\n", + i0,i1,i2,lamda01,lamda02,lamda12,v[0]); } /* ---------------------------------------------------------------------- @@ -2265,91 +2365,6 @@ int FixShake::unpack_exchange(int nlocal, double *buf) return m; } -/* ---------------------------------------------------------------------- - enforce SHAKE constraints from rRESPA - prediction portion is different than Verlet - rRESPA updating of atom coords is done with full v, but only portions of f -------------------------------------------------------------------------- */ - -void FixShake::post_force_respa(int vflag, int ilevel, int iloop) -{ - // call stats only on outermost level - - if (ilevel == nlevels_respa-1 && update->ntimestep == next_output) stats(); - - // perform SHAKE on every loop iteration of every rRESPA level - // except last loop iteration of inner levels - - if (ilevel < nlevels_respa-1 && iloop == loop_respa[ilevel]-1) return; - - // xshake = atom coords after next x update in innermost loop - // depends on rRESPA level - // for levels > 0 this includes more than one velocity update - // xshake = predicted position from call to this routine at level N = - // x + dt0 (v + dtN/m fN + 1/2 dt(N-1)/m f(N-1) + ... + 1/2 dt0/m f0) - - double ***f_level = ((FixRespa *) modify->fix[ifix_respa])->f_level; - dtfsq = dtf_inner * step_respa[ilevel]; - - double invmass,dtfmsq; - int jlevel; - - if (rmass) { - for (int i = 0; i < nlocal; i++) { - if (shake_flag[i]) { - invmass = 1.0 / rmass[i]; - dtfmsq = dtfsq * invmass; - xshake[i][0] = x[i][0] + dtv*v[i][0] + dtfmsq*f[i][0]; - xshake[i][1] = x[i][1] + dtv*v[i][1] + dtfmsq*f[i][1]; - xshake[i][2] = x[i][2] + dtv*v[i][2] + dtfmsq*f[i][2]; - for (jlevel = 0; jlevel < ilevel; jlevel++) { - dtfmsq = dtf_innerhalf * step_respa[jlevel] * invmass; - xshake[i][0] += dtfmsq*f_level[i][jlevel][0]; - xshake[i][1] += dtfmsq*f_level[i][jlevel][1]; - xshake[i][2] += dtfmsq*f_level[i][jlevel][2]; - } - } else xshake[i][2] = xshake[i][1] = xshake[i][0] = 0.0; - } - - } else { - for (int i = 0; i < nlocal; i++) { - if (shake_flag[i]) { - invmass = 1.0 / mass[type[i]]; - dtfmsq = dtfsq * invmass; - xshake[i][0] = x[i][0] + dtv*v[i][0] + dtfmsq*f[i][0]; - xshake[i][1] = x[i][1] + dtv*v[i][1] + dtfmsq*f[i][1]; - xshake[i][2] = x[i][2] + dtv*v[i][2] + dtfmsq*f[i][2]; - for (jlevel = 0; jlevel < ilevel; jlevel++) { - dtfmsq = dtf_innerhalf * step_respa[jlevel] * invmass; - xshake[i][0] += dtfmsq*f_level[i][jlevel][0]; - xshake[i][1] += dtfmsq*f_level[i][jlevel][1]; - xshake[i][2] += dtfmsq*f_level[i][jlevel][2]; - } - } else xshake[i][2] = xshake[i][1] = xshake[i][0] = 0.0; - } - } - - // communicate results if necessary - - if (nprocs > 1) comm->forward_comm_fix(this); - - // virial setup - - if (vflag) v_setup(vflag); - else evflag = 0; - - // loop over clusters - - int m; - for (int i = 0; i < nlist; i++) { - m = list[i]; - if (shake_flag[m] == 2) shake2(m); - else if (shake_flag[m] == 3) shake3(m); - else if (shake_flag[m] == 4) shake4(m); - else shake3angle(m); - } -} - /* ---------------------------------------------------------------------- */ int FixShake::pack_comm(int n, int *list, double *buf, int pbc_flag, int *pbc) diff --git a/src/fix_shake.h b/src/fix_shake.h index 8653cdd5cf..ebef049dbb 100644 --- a/src/fix_shake.h +++ b/src/fix_shake.h @@ -102,6 +102,7 @@ class FixShake : public Fix { void find_clusters(); int masscheck(double); void unconstrained_update(); + void unconstrained_update_respa(int); void shake2(int); void shake3(int); void shake4(int); diff --git a/src/lattice.h b/src/lattice.h index 358206c829..f38d2ac556 100644 --- a/src/lattice.h +++ b/src/lattice.h @@ -22,8 +22,8 @@ class Lattice : protected Pointers { public: int style; // enum list of NONE,SC,FCC,etc double xlattice,ylattice,zlattice; // lattice scale factors in 3 dims - double a1[3],a2[3],a3[3]; // vectors that bound unit cell - int nbasis; // # of atoms in basis of unit cell + double a1[3],a2[3],a3[3]; // edge vectors of unit cell + int nbasis; // # of basis atoms in unit cell double **basis; // fractional coords of each basis atom // within unit cell (0 <= coord < 1) diff --git a/src/modify.cpp b/src/modify.cpp index c0e9dc2e57..6da28a9772 100644 --- a/src/modify.cpp +++ b/src/modify.cpp @@ -348,21 +348,21 @@ double Modify::thermo_energy() 1st half of rRESPA integrate call, only for relevant fixes ------------------------------------------------------------------------- */ -void Modify::initial_integrate_respa(int vflag, int ilevel, int flag) +void Modify::initial_integrate_respa(int vflag, int ilevel, int iloop) { for (int i = 0; i < n_initial_integrate_respa; i++) fix[list_initial_integrate_respa[i]]-> - initial_integrate_respa(vflag,ilevel,flag); + initial_integrate_respa(vflag,ilevel,iloop); } /* ---------------------------------------------------------------------- rRESPA post_integrate call, only for relevant fixes ------------------------------------------------------------------------- */ -void Modify::post_integrate_respa(int ilevel, int flag) +void Modify::post_integrate_respa(int ilevel, int iloop) { for (int i = 0; i < n_post_integrate_respa; i++) - fix[list_post_integrate_respa[i]]->post_integrate_respa(ilevel,flag); + fix[list_post_integrate_respa[i]]->post_integrate_respa(ilevel,iloop); } /* ---------------------------------------------------------------------- @@ -389,10 +389,10 @@ void Modify::post_force_respa(int vflag, int ilevel, int iloop) 2nd half of rRESPA integrate call, only for relevant fixes ------------------------------------------------------------------------- */ -void Modify::final_integrate_respa(int ilevel) +void Modify::final_integrate_respa(int ilevel, int iloop) { for (int i = 0; i < n_final_integrate_respa; i++) - fix[list_final_integrate_respa[i]]->final_integrate_respa(ilevel); + fix[list_final_integrate_respa[i]]->final_integrate_respa(ilevel,iloop); } /* ---------------------------------------------------------------------- @@ -533,8 +533,10 @@ void Modify::add_fix(int narg, char **arg) // error if new style does not match old style // since can't replace it (all when-to-invoke ptrs would be invalid) // warn if new group != old group - // delete old fix - // set ptr to NULL in case new fix scans list of fixes + // delete old fix, but do not call update_callback(), + // since will replace this fix and thus other fix locs will not change + // set ptr to NULL in case new fix scans list of fixes, + // e.g. scan will occur in add_callback() if called by new fix // if fix ID does not exist: // set newflag = 1 so create new fix // extend fix and fmask lists as necessary diff --git a/src/modify.h b/src/modify.h index 07d53116e9..03ad635996 100644 --- a/src/modify.h +++ b/src/modify.h @@ -53,11 +53,11 @@ class Modify : protected Pointers { void end_of_step(); double thermo_energy(); - void initial_integrate_respa(int,int,int); - void post_integrate_respa(int,int); - void pre_force_respa(int,int,int); - void post_force_respa(int,int,int); - void final_integrate_respa(int); + void initial_integrate_respa(int, int, int); + void post_integrate_respa(int, int); + void pre_force_respa(int, int, int); + void post_force_respa(int, int, int); + void final_integrate_respa(int, int); void min_pre_exchange(); void min_post_force(int); diff --git a/src/respa.cpp b/src/respa.cpp index 8d89cbbe13..f898d88d7d 100644 --- a/src/respa.cpp +++ b/src/respa.cpp @@ -475,24 +475,17 @@ void Respa::recurse(int ilevel) for (int iloop = 0; iloop < loop[ilevel]; iloop++) { - modify->initial_integrate_respa(vflag,ilevel,0); + modify->initial_integrate_respa(vflag,ilevel,iloop); if (modify->n_post_integrate_respa) modify->post_integrate_respa(ilevel,iloop); if (ilevel) recurse(ilevel-1); - // at outermost level: - // call initial_integrate w/ flag = 1 so fix NPT,NPH performs - // 2nd box deform at correct symmetric position in rRESPA hierarchy - // check on rebuilding neighbor list + // at outermost level, check on rebuilding neighbor list // at innermost level, communicate // at middle levels, do nothing if (ilevel == nlevels-1) { - modify->initial_integrate_respa(vflag,ilevel,1); - if (modify->n_post_integrate_respa) - modify->post_integrate_respa(ilevel,iloop); - int nflag = neighbor->decide(); if (nflag) { if (modify->n_pre_exchange) modify->pre_exchange(); @@ -570,7 +563,7 @@ void Respa::recurse(int ilevel) if (modify->n_post_force_respa) modify->post_force_respa(vflag,ilevel,iloop); - modify->final_integrate_respa(ilevel); + modify->final_integrate_respa(ilevel,iloop); } copy_f_flevel(ilevel);