diff --git a/src/ASPHERE/fix_npt_asphere.cpp b/src/ASPHERE/fix_npt_asphere.cpp index 6fd81afe92..4c1403f83d 100755 --- a/src/ASPHERE/fix_npt_asphere.cpp +++ b/src/ASPHERE/fix_npt_asphere.cpp @@ -43,7 +43,7 @@ FixNPTASphere::FixNPTASphere(LAMMPS *lmp, int narg, char **arg) : 1st half of Verlet update ------------------------------------------------------------------------- */ -void FixNPTASphere::initial_integrate() +void FixNPTASphere::initial_integrate(int vflag) { int i; diff --git a/src/ASPHERE/fix_npt_asphere.h b/src/ASPHERE/fix_npt_asphere.h index beec20f923..8cad21b0bc 100755 --- a/src/ASPHERE/fix_npt_asphere.h +++ b/src/ASPHERE/fix_npt_asphere.h @@ -22,7 +22,7 @@ class FixNPTASphere : public FixNPT { public: FixNPTASphere(class LAMMPS *, int, char **); ~FixNPTASphere() {} - void initial_integrate(); + void initial_integrate(int); void final_integrate(); private: diff --git a/src/ASPHERE/fix_nve_asphere.cpp b/src/ASPHERE/fix_nve_asphere.cpp index 3972967afc..722a8d0838 100755 --- a/src/ASPHERE/fix_nve_asphere.cpp +++ b/src/ASPHERE/fix_nve_asphere.cpp @@ -60,7 +60,7 @@ void FixNVEASphere::init() /* ---------------------------------------------------------------------- */ -void FixNVEASphere::initial_integrate() +void FixNVEASphere::initial_integrate(int vflag) { double dtfm; diff --git a/src/ASPHERE/fix_nve_asphere.h b/src/ASPHERE/fix_nve_asphere.h index b41a4c0e51..c748b6db5e 100755 --- a/src/ASPHERE/fix_nve_asphere.h +++ b/src/ASPHERE/fix_nve_asphere.h @@ -23,7 +23,7 @@ class FixNVEASphere : public FixNVE { FixNVEASphere(class LAMMPS *, int, char **); ~FixNVEASphere(); void init(); - void initial_integrate(); + void initial_integrate(int); void final_integrate(); private: diff --git a/src/ASPHERE/fix_nvt_asphere.cpp b/src/ASPHERE/fix_nvt_asphere.cpp index fad60a7e73..5fec36e1fe 100755 --- a/src/ASPHERE/fix_nvt_asphere.cpp +++ b/src/ASPHERE/fix_nvt_asphere.cpp @@ -43,7 +43,7 @@ FixNVTASphere::FixNVTASphere(LAMMPS *lmp, int narg, char **arg) : /* ---------------------------------------------------------------------- */ -void FixNVTASphere::initial_integrate() +void FixNVTASphere::initial_integrate(int vflag) { double dtfm; diff --git a/src/ASPHERE/fix_nvt_asphere.h b/src/ASPHERE/fix_nvt_asphere.h index 6a12c65dd3..8f037d541b 100755 --- a/src/ASPHERE/fix_nvt_asphere.h +++ b/src/ASPHERE/fix_nvt_asphere.h @@ -22,7 +22,7 @@ class FixNVTASphere : public FixNVT { public: FixNVTASphere(class LAMMPS *, int, char **); ~FixNVTASphere() {} - void initial_integrate(); + void initial_integrate(int); void final_integrate(); private: diff --git a/src/DIPOLE/fix_nve_dipole.cpp b/src/DIPOLE/fix_nve_dipole.cpp index 2048e3bbe2..8cbccf9541 100644 --- a/src/DIPOLE/fix_nve_dipole.cpp +++ b/src/DIPOLE/fix_nve_dipole.cpp @@ -79,7 +79,7 @@ void FixNVEDipole::init() /* ---------------------------------------------------------------------- */ -void FixNVEDipole::initial_integrate() +void FixNVEDipole::initial_integrate(int vflag) { double dtfm,msq,scale; @@ -169,14 +169,14 @@ void FixNVEDipole::final_integrate() /* ---------------------------------------------------------------------- */ -void FixNVEDipole::initial_integrate_respa(int ilevel, int flag) +void FixNVEDipole::initial_integrate_respa(int vflag, int ilevel, int flag) { if (flag) return; // only used by NPT,NPH dtv = step_respa[ilevel]; dtf = 0.5 * step_respa[ilevel] * force->ftm2v; - if (ilevel == 0) initial_integrate(); + if (ilevel == 0) initial_integrate(vflag); else final_integrate(); } diff --git a/src/DIPOLE/fix_nve_dipole.h b/src/DIPOLE/fix_nve_dipole.h index 01b3f75fda..3c3417832a 100644 --- a/src/DIPOLE/fix_nve_dipole.h +++ b/src/DIPOLE/fix_nve_dipole.h @@ -24,9 +24,9 @@ class FixNVEDipole : public Fix { ~FixNVEDipole(); int setmask(); void init(); - void initial_integrate(); + void initial_integrate(int); void final_integrate(); - void initial_integrate_respa(int,int); + void initial_integrate_respa(int, int, int); void final_integrate_respa(int); void reset_dt(); diff --git a/src/GRANULAR/fix_freeze.cpp b/src/GRANULAR/fix_freeze.cpp index d1f7e4a0f8..67e3d1dbd8 100644 --- a/src/GRANULAR/fix_freeze.cpp +++ b/src/GRANULAR/fix_freeze.cpp @@ -55,9 +55,9 @@ void FixFreeze::init() /* ---------------------------------------------------------------------- */ -void FixFreeze::setup() +void FixFreeze::setup(int vflag) { - post_force(1); + post_force(vflag); } /* ---------------------------------------------------------------------- */ diff --git a/src/GRANULAR/fix_freeze.h b/src/GRANULAR/fix_freeze.h index e23d3ca005..a003e79b62 100644 --- a/src/GRANULAR/fix_freeze.h +++ b/src/GRANULAR/fix_freeze.h @@ -23,7 +23,7 @@ class FixFreeze : public Fix { FixFreeze(class LAMMPS *, int, char **); int setmask(); void init(); - void setup(); + void setup(int); void post_force(int); }; diff --git a/src/GRANULAR/fix_nve_gran.cpp b/src/GRANULAR/fix_nve_gran.cpp index 1afa347071..7b52b633e6 100644 --- a/src/GRANULAR/fix_nve_gran.cpp +++ b/src/GRANULAR/fix_nve_gran.cpp @@ -61,7 +61,7 @@ void FixNVEGran::init() /* ---------------------------------------------------------------------- */ -void FixNVEGran::initial_integrate() +void FixNVEGran::initial_integrate(int vflag) { double dtfm; diff --git a/src/GRANULAR/fix_nve_gran.h b/src/GRANULAR/fix_nve_gran.h index ebe1e73828..4aa111ffb8 100644 --- a/src/GRANULAR/fix_nve_gran.h +++ b/src/GRANULAR/fix_nve_gran.h @@ -23,7 +23,7 @@ class FixNVEGran : public Fix { FixNVEGran(class LAMMPS *, int, char **); int setmask(); void init(); - void initial_integrate(); + void initial_integrate(int); void final_integrate(); void reset_dt(); diff --git a/src/GRANULAR/fix_wall_gran.cpp b/src/GRANULAR/fix_wall_gran.cpp index 5b123becb9..64a8e970c6 100644 --- a/src/GRANULAR/fix_wall_gran.cpp +++ b/src/GRANULAR/fix_wall_gran.cpp @@ -196,9 +196,9 @@ void FixWallGran::init() /* ---------------------------------------------------------------------- */ -void FixWallGran::setup() +void FixWallGran::setup(int vflag) { - post_force(1); + post_force(vflag); } /* ---------------------------------------------------------------------- */ diff --git a/src/GRANULAR/fix_wall_gran.h b/src/GRANULAR/fix_wall_gran.h index 0fe97a7dfa..cae688cd13 100644 --- a/src/GRANULAR/fix_wall_gran.h +++ b/src/GRANULAR/fix_wall_gran.h @@ -24,7 +24,7 @@ class FixWallGran : public Fix { ~FixWallGran(); int setmask(); void init(); - void setup(); + void setup(int); void post_force(int); double memory_usage(); diff --git a/src/compute_stress_atom.cpp b/src/compute_stress_atom.cpp index 86d89897f9..0e472ab90f 100644 --- a/src/compute_stress_atom.cpp +++ b/src/compute_stress_atom.cpp @@ -22,6 +22,8 @@ #include "angle.h" #include "dihedral.h" #include "improper.h" +#include "modify.h" +#include "fix.h" #include "memory.h" #include "error.h" @@ -46,10 +48,12 @@ ComputeStressAtom::ComputeStressAtom(LAMMPS *lmp, int narg, char **arg) : keflag = 1; pairflag = 1; bondflag = angleflag = dihedralflag = improperflag = 1; + fixflag = 1; } else { keflag = 0; pairflag = 0; bondflag = angleflag = dihedralflag = improperflag = 0; + fixflag = 0; int iarg = 3; while (iarg < narg) { if (strcmp(arg[iarg],"ke") == 0) keflag = 1; @@ -58,6 +62,7 @@ ComputeStressAtom::ComputeStressAtom(LAMMPS *lmp, int narg, char **arg) : else if (strcmp(arg[iarg],"angle") == 0) angleflag = 1; else if (strcmp(arg[iarg],"dihedral") == 0) dihedralflag = 1; else if (strcmp(arg[iarg],"improper") == 0) improperflag = 1; + else if (strcmp(arg[iarg],"fix") == 0) fixflag = 1; else error->all("Illegal compute stress/atom command"); iarg++; } @@ -148,7 +153,19 @@ void ComputeStressAtom::compute_peratom() stress[i][j] += vatom[i][j]; } - // communicate ghost energy between neighbor procs + // add in per-atom contributions from relevant fixes + + if (fixflag) { + for (int i = 0; i < modify->nfix; i++) + if (modify->fix[i]->virial_flag) { + double **vatom = modify->fix[i]->vatom; + for (i = 0; i < nlocal; i++) + for (j = 0; j < 6; j++) + stress[i][j] += vatom[i][j]; + } + } + + // communicate ghost atom virials between neighbor procs if (force->newton) comm->reverse_comm_compute(this); diff --git a/src/compute_stress_atom.h b/src/compute_stress_atom.h index d1cff5530b..98b727bf80 100644 --- a/src/compute_stress_atom.h +++ b/src/compute_stress_atom.h @@ -29,7 +29,7 @@ class ComputeStressAtom : public Compute { double memory_usage(); private: - int keflag,pairflag,bondflag,angleflag,dihedralflag,improperflag; + int keflag,pairflag,bondflag,angleflag,dihedralflag,improperflag,fixflag; int nmax; double **stress; }; diff --git a/src/fix.cpp b/src/fix.cpp index dafb55c21c..b822fdf1fb 100644 --- a/src/fix.cpp +++ b/src/fix.cpp @@ -14,7 +14,9 @@ #include "string.h" #include "ctype.h" #include "fix.h" +#include "atom.h" #include "group.h" +#include "memory.h" #include "error.h" using namespace LAMMPS_NS; @@ -55,6 +57,9 @@ Fix::Fix(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) comm_forward = comm_reverse = 0; + maxvatom = 0; + vatom = NULL; + // mask settings - same as in modify.cpp INITIAL_INTEGRATE = 1; @@ -77,6 +82,7 @@ Fix::~Fix() { delete [] id; delete [] style; + memory->destroy_2d_double_array(vatom); } /* ---------------------------------------------------------------------- @@ -103,3 +109,78 @@ void Fix::modify_params(int narg, char **arg) } } } + +/* ---------------------------------------------------------------------- + setup for virial computation + see integrate::ev_set() for values of vflag (0-6) +------------------------------------------------------------------------- */ + +void Fix::v_setup(int vflag) +{ + int i,n; + + evflag = 1; + + vflag_global = vflag % 4; + vflag_atom = vflag / 4; + + // reallocate per-atom array if necessary + + if (vflag_atom && atom->nlocal > maxvatom) { + maxvatom = atom->nmax; + memory->destroy_2d_double_array(vatom); + vatom = memory->create_2d_double_array(maxvatom,6,"bond:vatom"); + } + + // zero accumulators + + if (vflag_global) for (i = 0; i < 6; i++) virial[i] = 0.0; + if (vflag_atom) { + n = atom->nlocal; + for (i = 0; i < n; i++) { + vatom[i][0] = 0.0; + vatom[i][1] = 0.0; + vatom[i][2] = 0.0; + vatom[i][3] = 0.0; + vatom[i][4] = 0.0; + vatom[i][5] = 0.0; + } + } +} + +/* ---------------------------------------------------------------------- + tally virial into global and per-atom accumulators + v = total virial for the interaction involving total atoms + n = # of local atoms involved, with local indices in list + increment global virial by n/total fraction + increment per-atom virial of each atom in list by 1/total fraction + assumes other procs will tally left-over fractions +------------------------------------------------------------------------- */ + +void Fix::v_tally(int n, int *list, double total, double *v) +{ + int m; + + if (vflag_global) { + double fraction = n/total; + virial[0] += fraction*v[0]; + virial[1] += fraction*v[1]; + virial[2] += fraction*v[2]; + virial[3] += fraction*v[3]; + virial[4] += fraction*v[4]; + virial[5] += fraction*v[5]; + } + + if (vflag_atom) { + double fraction = 1.0/total; + for (int i = 0; i < n; i++) { + m = list[i]; + vatom[m][0] += fraction*v[0]; + vatom[m][1] += fraction*v[1]; + vatom[m][2] += fraction*v[2]; + vatom[m][3] += fraction*v[3]; + vatom[m][4] += fraction*v[4]; + vatom[m][5] += fraction*v[5]; + } + } +} diff --git a/src/fix.h b/src/fix.h index 754daf3c4c..a83886639a 100644 --- a/src/fix.h +++ b/src/fix.h @@ -51,7 +51,8 @@ class Fix : protected Pointers { int comm_forward; // size of forward communication (0 if none) int comm_reverse; // size of reverse communication (0 if none) - double virial[6]; // fix contribution to pressure virial + double virial[6]; // accumlated virial + double **vatom; // accumulated per-atom virial int INITIAL_INTEGRATE,PRE_EXCHANGE,PRE_NEIGHBOR; // mask settings int POST_FORCE,FINAL_INTEGRATE,END_OF_STEP,THERMO_ENERGY; @@ -66,9 +67,9 @@ class Fix : protected Pointers { virtual void init() {} virtual void init_list(int, class NeighList *) {} - virtual void setup() {} - virtual void min_setup() {} - virtual void initial_integrate() {} + virtual void setup(int) {} + virtual void min_setup(int) {} + virtual void initial_integrate(int) {} virtual void pre_exchange() {} virtual void pre_neighbor() {} virtual void post_force(int) {} @@ -86,7 +87,7 @@ class Fix : protected Pointers { virtual int size_restart(int) {return 0;} virtual int maxsize_restart() {return 0;} - virtual void initial_integrate_respa(int, int) {} + virtual void initial_integrate_respa(int, int, int) {} virtual void post_force_respa(int, int, int) {} virtual void final_integrate_respa(int) {} @@ -110,6 +111,14 @@ class Fix : protected Pointers { virtual int modify_param(int, char **) {return 0;} virtual double memory_usage() {return 0.0;} + + protected: + int evflag; + int vflag_global,vflag_atom; + int maxvatom; + + void v_setup(int); + void v_tally(int, int *, double, double *); }; } diff --git a/src/fix_add_force.cpp b/src/fix_add_force.cpp index a865b9a63a..f5bef046e4 100644 --- a/src/fix_add_force.cpp +++ b/src/fix_add_force.cpp @@ -53,22 +53,22 @@ void FixAddForce::init() /* ---------------------------------------------------------------------- */ -void FixAddForce::setup() +void FixAddForce::setup(int vflag) { if (strcmp(update->integrate_style,"verlet") == 0) - post_force(1); + post_force(vflag); else { ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1); - post_force_respa(1,nlevels_respa-1,0); + post_force_respa(vflag,nlevels_respa-1,0); ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1); } } /* ---------------------------------------------------------------------- */ -void FixAddForce::min_setup() +void FixAddForce::min_setup(int vflag) { - post_force(1); + post_force(vflag); } /* ---------------------------------------------------------------------- */ diff --git a/src/fix_add_force.h b/src/fix_add_force.h index 44d68de343..0253ad0319 100644 --- a/src/fix_add_force.h +++ b/src/fix_add_force.h @@ -23,8 +23,8 @@ class FixAddForce : public Fix { FixAddForce(class LAMMPS *, int, char **); int setmask(); void init(); - void setup(); - void min_setup(); + void setup(int); + void min_setup(int); void post_force(int); void post_force_respa(int, int, int); void min_post_force(int); diff --git a/src/fix_ave_atom.cpp b/src/fix_ave_atom.cpp index bb320105c8..a5a1389f8e 100644 --- a/src/fix_ave_atom.cpp +++ b/src/fix_ave_atom.cpp @@ -255,7 +255,7 @@ void FixAveAtom::init() only does something if nvalid = current timestep ------------------------------------------------------------------------- */ -void FixAveAtom::setup() +void FixAveAtom::setup(int vflag) { end_of_step(); } diff --git a/src/fix_ave_atom.h b/src/fix_ave_atom.h index e2f47a417c..e54e0512d8 100644 --- a/src/fix_ave_atom.h +++ b/src/fix_ave_atom.h @@ -25,7 +25,7 @@ class FixAveAtom : public Fix { ~FixAveAtom(); int setmask(); void init(); - void setup(); + void setup(int); void end_of_step(); double memory_usage(); diff --git a/src/fix_ave_force.cpp b/src/fix_ave_force.cpp index e05393918e..318ca411a0 100644 --- a/src/fix_ave_force.cpp +++ b/src/fix_ave_force.cpp @@ -69,23 +69,23 @@ void FixAveForce::init() /* ---------------------------------------------------------------------- */ -void FixAveForce::setup() +void FixAveForce::setup(int vflag) { if (strcmp(update->integrate_style,"verlet") == 0) - post_force(1); + post_force(vflag); else for (int ilevel = 0; ilevel < nlevels_respa; ilevel++) { ((Respa *) update->integrate)->copy_flevel_f(ilevel); - post_force_respa(1,ilevel,0); + post_force_respa(vflag,ilevel,0); ((Respa *) update->integrate)->copy_f_flevel(ilevel); } } /* ---------------------------------------------------------------------- */ -void FixAveForce::min_setup() +void FixAveForce::min_setup(int vflag) { - post_force(1); + post_force(vflag); } /* ---------------------------------------------------------------------- */ diff --git a/src/fix_ave_force.h b/src/fix_ave_force.h index 542b49befa..b5b32d452f 100644 --- a/src/fix_ave_force.h +++ b/src/fix_ave_force.h @@ -23,8 +23,8 @@ class FixAveForce : public Fix { FixAveForce(class LAMMPS *, int, char **); int setmask(); void init(); - void setup(); - void min_setup(); + void setup(int); + void min_setup(int); void post_force(int); void post_force_respa(int, int, int); void min_post_force(int); diff --git a/src/fix_ave_spatial.cpp b/src/fix_ave_spatial.cpp index 093e6ed573..bcdac37402 100644 --- a/src/fix_ave_spatial.cpp +++ b/src/fix_ave_spatial.cpp @@ -403,7 +403,7 @@ void FixAveSpatial::init() only does something if nvalid = current timestep ------------------------------------------------------------------------- */ -void FixAveSpatial::setup() +void FixAveSpatial::setup(int vflag) { end_of_step(); } diff --git a/src/fix_ave_spatial.h b/src/fix_ave_spatial.h index fd2b089116..0a4da9876e 100644 --- a/src/fix_ave_spatial.h +++ b/src/fix_ave_spatial.h @@ -25,7 +25,7 @@ class FixAveSpatial : public Fix { ~FixAveSpatial(); int setmask(); void init(); - void setup(); + void setup(int); void end_of_step(); double compute_vector(int); double memory_usage(); diff --git a/src/fix_ave_time.cpp b/src/fix_ave_time.cpp index a0bc80df02..f2fd1936dc 100644 --- a/src/fix_ave_time.cpp +++ b/src/fix_ave_time.cpp @@ -314,7 +314,7 @@ void FixAveTime::init() only does something if nvalid = current timestep ------------------------------------------------------------------------- */ -void FixAveTime::setup() +void FixAveTime::setup(int vflag) { end_of_step(); } diff --git a/src/fix_ave_time.h b/src/fix_ave_time.h index ea0103e9ea..7311a72fcf 100644 --- a/src/fix_ave_time.h +++ b/src/fix_ave_time.h @@ -25,7 +25,7 @@ class FixAveTime : public Fix { ~FixAveTime(); int setmask(); void init(); - void setup(); + void setup(int); void end_of_step(); double compute_scalar(); double compute_vector(int); diff --git a/src/fix_com.cpp b/src/fix_com.cpp index c52f78db5a..689d48505b 100644 --- a/src/fix_com.cpp +++ b/src/fix_com.cpp @@ -72,7 +72,7 @@ void FixCOM::init() /* ---------------------------------------------------------------------- */ -void FixCOM::setup() +void FixCOM::setup(int vflag) { if (first) end_of_step(); first = 0; diff --git a/src/fix_com.h b/src/fix_com.h index 4623e0471a..f936508b17 100644 --- a/src/fix_com.h +++ b/src/fix_com.h @@ -25,7 +25,7 @@ class FixCOM : public Fix { ~FixCOM(); int setmask(); void init(); - void setup(); + void setup(int); void end_of_step(); private: diff --git a/src/fix_drag.cpp b/src/fix_drag.cpp index 72fdf9c8db..0d82374340 100644 --- a/src/fix_drag.cpp +++ b/src/fix_drag.cpp @@ -71,13 +71,13 @@ void FixDrag::init() /* ---------------------------------------------------------------------- */ -void FixDrag::setup() +void FixDrag::setup(int vflag) { if (strcmp(update->integrate_style,"verlet") == 0) - post_force(1); + post_force(vflag); else { ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1); - post_force_respa(1,nlevels_respa-1,0); + post_force_respa(vflag,nlevels_respa-1,0); ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1); } } diff --git a/src/fix_drag.h b/src/fix_drag.h index 65c074f73c..6db9ba1114 100644 --- a/src/fix_drag.h +++ b/src/fix_drag.h @@ -23,7 +23,7 @@ class FixDrag : public Fix { FixDrag(class LAMMPS *, int, char **); int setmask(); void init(); - void setup(); + void setup(int); void post_force(int); void post_force_respa(int, int, int); double compute_vector(int); diff --git a/src/fix_dt_reset.cpp b/src/fix_dt_reset.cpp index bab0d0e6b7..358e770699 100644 --- a/src/fix_dt_reset.cpp +++ b/src/fix_dt_reset.cpp @@ -120,7 +120,7 @@ void FixDtReset::init() /* ---------------------------------------------------------------------- */ -void FixDtReset::setup() +void FixDtReset::setup(int vflag) { end_of_step(); } diff --git a/src/fix_dt_reset.h b/src/fix_dt_reset.h index fa0369dc44..7099c33a3b 100644 --- a/src/fix_dt_reset.h +++ b/src/fix_dt_reset.h @@ -24,7 +24,7 @@ class FixDtReset : public Fix { ~FixDtReset() {} int setmask(); void init(); - void setup(); + void setup(int); void end_of_step(); double compute_scalar(); double compute_vector(int); diff --git a/src/fix_efield.cpp b/src/fix_efield.cpp index 96d3029d97..a3fb43eee2 100644 --- a/src/fix_efield.cpp +++ b/src/fix_efield.cpp @@ -65,13 +65,13 @@ void FixEfield::init() /* ---------------------------------------------------------------------- */ -void FixEfield::setup() +void FixEfield::setup(int vflag) { if (strcmp(update->integrate_style,"verlet") == 0) - post_force(1); + post_force(vflag); else { ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1); - post_force_respa(1,nlevels_respa-1,0); + post_force_respa(vflag,nlevels_respa-1,0); ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1); } } diff --git a/src/fix_efield.h b/src/fix_efield.h index ce973abab6..3d42822a04 100644 --- a/src/fix_efield.h +++ b/src/fix_efield.h @@ -23,7 +23,7 @@ class FixEfield : public Fix { FixEfield(class LAMMPS *, int, char **); int setmask(); void init(); - void setup(); + void setup(int); void post_force(int); void post_force_respa(int, int, int); diff --git a/src/fix_enforce2d.cpp b/src/fix_enforce2d.cpp index 2ca68ed17e..fff0f53535 100644 --- a/src/fix_enforce2d.cpp +++ b/src/fix_enforce2d.cpp @@ -41,15 +41,15 @@ int FixEnforce2D::setmask() /* ---------------------------------------------------------------------- */ -void FixEnforce2D::setup() +void FixEnforce2D::setup(int vflag) { if (strcmp(update->integrate_style,"verlet") == 0) - post_force(1); + post_force(vflag); else { int nlevels_respa = ((Respa *) update->integrate)->nlevels; for (int ilevel = 0; ilevel < nlevels_respa; ilevel++) { ((Respa *) update->integrate)->copy_flevel_f(ilevel); - post_force_respa(1,ilevel,0); + post_force_respa(vflag,ilevel,0); ((Respa *) update->integrate)->copy_f_flevel(ilevel); } } @@ -57,9 +57,9 @@ void FixEnforce2D::setup() /* ---------------------------------------------------------------------- */ -void FixEnforce2D::min_setup() +void FixEnforce2D::min_setup(int vflag) { - post_force(1); + post_force(vflag); } /* ---------------------------------------------------------------------- */ diff --git a/src/fix_enforce2d.h b/src/fix_enforce2d.h index eade63b0da..b1b9440e72 100644 --- a/src/fix_enforce2d.h +++ b/src/fix_enforce2d.h @@ -22,8 +22,8 @@ class FixEnforce2D : public Fix { public: FixEnforce2D(class LAMMPS *, int, char **); int setmask(); - void setup(); - void min_setup(); + void setup(int); + void min_setup(int); void post_force(int); void post_force_respa(int, int, int); void min_post_force(int); diff --git a/src/fix_gravity.cpp b/src/fix_gravity.cpp index 90ac0247b8..576b1d14b3 100644 --- a/src/fix_gravity.cpp +++ b/src/fix_gravity.cpp @@ -111,9 +111,9 @@ void FixGravity::init() /* ---------------------------------------------------------------------- */ -void FixGravity::setup() +void FixGravity::setup(int vflag) { - post_force(1); + post_force(vflag); } /* ---------------------------------------------------------------------- */ diff --git a/src/fix_gravity.h b/src/fix_gravity.h index 027595d1ad..017a5246ff 100644 --- a/src/fix_gravity.h +++ b/src/fix_gravity.h @@ -25,7 +25,7 @@ class FixGravity : public Fix { FixGravity(class LAMMPS *, int, char **); int setmask(); void init(); - void setup(); + void setup(int); void post_force(int); private: diff --git a/src/fix_gyration.cpp b/src/fix_gyration.cpp index d21a91fd8e..0b8ed1668a 100644 --- a/src/fix_gyration.cpp +++ b/src/fix_gyration.cpp @@ -72,7 +72,7 @@ void FixGyration::init() /* ---------------------------------------------------------------------- */ -void FixGyration::setup() +void FixGyration::setup(int vflag) { if (first) end_of_step(); first = 0; diff --git a/src/fix_gyration.h b/src/fix_gyration.h index 448a80cd06..b0b52450ee 100644 --- a/src/fix_gyration.h +++ b/src/fix_gyration.h @@ -25,7 +25,7 @@ class FixGyration : public Fix { ~FixGyration(); int setmask(); void init(); - void setup(); + void setup(int); void end_of_step(); private: diff --git a/src/fix_indent.cpp b/src/fix_indent.cpp index 9c697e7546..72a1ef47ab 100644 --- a/src/fix_indent.cpp +++ b/src/fix_indent.cpp @@ -131,22 +131,22 @@ void FixIndent::init() /* ---------------------------------------------------------------------- */ -void FixIndent::setup() +void FixIndent::setup(int vflag) { if (strcmp(update->integrate_style,"verlet") == 0) - post_force(1); + post_force(vflag); else { ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1); - post_force_respa(1,nlevels_respa-1,0); + post_force_respa(vflag,nlevels_respa-1,0); ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1); } } /* ---------------------------------------------------------------------- */ -void FixIndent::min_setup() +void FixIndent::min_setup(int vflag) { - post_force(1); + post_force(vflag); } /* ---------------------------------------------------------------------- */ diff --git a/src/fix_indent.h b/src/fix_indent.h index f0604146da..ce092666b1 100644 --- a/src/fix_indent.h +++ b/src/fix_indent.h @@ -23,8 +23,8 @@ class FixIndent : public Fix { FixIndent(class LAMMPS *, int, char **); int setmask(); void init(); - void setup(); - void min_setup(); + void setup(int); + void min_setup(int); virtual void post_force(int); void post_force_respa(int, int, int); void min_post_force(int); diff --git a/src/fix_langevin.cpp b/src/fix_langevin.cpp index 48f7a3c9ea..9a31cb89fd 100644 --- a/src/fix_langevin.cpp +++ b/src/fix_langevin.cpp @@ -126,13 +126,13 @@ void FixLangevin::init() /* ---------------------------------------------------------------------- */ -void FixLangevin::setup() +void FixLangevin::setup(int vflag) { if (strcmp(update->integrate_style,"verlet") == 0) - post_force(1); + post_force(vflag); else { ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1); - post_force_respa(1,nlevels_respa-1,0); + post_force_respa(vflag,nlevels_respa-1,0); ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1); } } diff --git a/src/fix_langevin.h b/src/fix_langevin.h index 8a8a89aeb8..1deeb84926 100644 --- a/src/fix_langevin.h +++ b/src/fix_langevin.h @@ -24,7 +24,7 @@ class FixLangevin : public Fix { ~FixLangevin(); int setmask(); void init(); - void setup(); + void setup(int); void post_force(int); void post_force_respa(int, int, int); void reset_target(double); diff --git a/src/fix_line_force.cpp b/src/fix_line_force.cpp index bde78c5040..d6750b4d01 100644 --- a/src/fix_line_force.cpp +++ b/src/fix_line_force.cpp @@ -45,15 +45,15 @@ int FixLineForce::setmask() /* ---------------------------------------------------------------------- */ -void FixLineForce::setup() +void FixLineForce::setup(int vflag) { if (strcmp(update->integrate_style,"verlet") == 0) - post_force(1); + post_force(vflag); else { int nlevels_respa = ((Respa *) update->integrate)->nlevels; for (int ilevel = 0; ilevel < nlevels_respa; ilevel++) { ((Respa *) update->integrate)->copy_flevel_f(ilevel); - post_force_respa(1,ilevel,0); + post_force_respa(vflag,ilevel,0); ((Respa *) update->integrate)->copy_f_flevel(ilevel); } } @@ -61,9 +61,9 @@ void FixLineForce::setup() /* ---------------------------------------------------------------------- */ -void FixLineForce::min_setup() +void FixLineForce::min_setup(int vflag) { - post_force(1); + post_force(vflag); } /* ---------------------------------------------------------------------- */ diff --git a/src/fix_line_force.h b/src/fix_line_force.h index e7fbc20b24..eed65e7458 100644 --- a/src/fix_line_force.h +++ b/src/fix_line_force.h @@ -22,8 +22,8 @@ class FixLineForce : public Fix { public: FixLineForce(class LAMMPS *, int, char **); int setmask(); - void setup(); - void min_setup(); + void setup(int); + void min_setup(int); void post_force(int); void post_force_respa(int, int, int); void min_post_force(int); diff --git a/src/fix_msd.cpp b/src/fix_msd.cpp index eb816de80e..f593c25cc8 100644 --- a/src/fix_msd.cpp +++ b/src/fix_msd.cpp @@ -123,7 +123,7 @@ void FixMSD::init() /* ---------------------------------------------------------------------- */ -void FixMSD::setup() +void FixMSD::setup(int vflag) { if (first) end_of_step(); first = 0; diff --git a/src/fix_msd.h b/src/fix_msd.h index edade856d8..07e1817e93 100644 --- a/src/fix_msd.h +++ b/src/fix_msd.h @@ -25,7 +25,7 @@ class FixMSD : public Fix { ~FixMSD(); int setmask(); void init(); - void setup(); + void setup(int); void end_of_step(); double memory_usage(); diff --git a/src/fix_nph.cpp b/src/fix_nph.cpp index 62dfe86138..0556f593d4 100644 --- a/src/fix_nph.cpp +++ b/src/fix_nph.cpp @@ -306,7 +306,7 @@ void FixNPH::init() compute T,P before integrator starts ------------------------------------------------------------------------- */ -void FixNPH::setup() +void FixNPH::setup(int vflag) { p_target[0] = p_start[0]; // used by compute_scalar() p_target[1] = p_start[1]; @@ -330,7 +330,7 @@ void FixNPH::setup() 1st half of Verlet update ------------------------------------------------------------------------- */ -void FixNPH::initial_integrate() +void FixNPH::initial_integrate(int vflag) { int i; @@ -455,7 +455,7 @@ void FixNPH::final_integrate() /* ---------------------------------------------------------------------- */ -void FixNPH::initial_integrate_respa(int ilevel, int flag) +void FixNPH::initial_integrate_respa(int vflag, int ilevel, int flag) { // if flag = 1, then is 2nd call at outermost level from rRESPA // perform 2nd half of box remap on own + ghost atoms and return diff --git a/src/fix_nph.h b/src/fix_nph.h index fe003a2d88..8560a6444e 100644 --- a/src/fix_nph.h +++ b/src/fix_nph.h @@ -24,10 +24,10 @@ class FixNPH : public Fix { ~FixNPH(); int setmask(); void init(); - void setup(); - void initial_integrate(); + void setup(int); + void initial_integrate(int); void final_integrate(); - void initial_integrate_respa(int,int); + void initial_integrate_respa(int, int, int); void final_integrate_respa(int); double compute_scalar(); void write_restart(FILE *); diff --git a/src/fix_npt.cpp b/src/fix_npt.cpp index 820ad96e36..7b51897139 100644 --- a/src/fix_npt.cpp +++ b/src/fix_npt.cpp @@ -307,7 +307,7 @@ void FixNPT::init() compute T,P before integrator starts ------------------------------------------------------------------------- */ -void FixNPT::setup() +void FixNPT::setup(int vflag) { t_target = t_start; // used by compute_scalar() p_target[0] = p_start[0]; @@ -332,7 +332,7 @@ void FixNPT::setup() 1st half of Verlet update ------------------------------------------------------------------------- */ -void FixNPT::initial_integrate() +void FixNPT::initial_integrate(int vflag) { int i; @@ -471,7 +471,7 @@ void FixNPT::final_integrate() /* ---------------------------------------------------------------------- */ -void FixNPT::initial_integrate_respa(int ilevel, int flag) +void FixNPT::initial_integrate_respa(int vflag, int ilevel, int flag) { // if flag = 1, then is 2nd call at outermost level from rRESPA // perform 2nd half of box remap on own + ghost atoms and return diff --git a/src/fix_npt.h b/src/fix_npt.h index 9505e7513e..f40370711f 100644 --- a/src/fix_npt.h +++ b/src/fix_npt.h @@ -24,10 +24,10 @@ class FixNPT : public Fix { virtual ~FixNPT(); int setmask(); void init(); - void setup(); - virtual void initial_integrate(); + void setup(int); + void initial_integrate(int); virtual void final_integrate(); - void initial_integrate_respa(int, int); + void initial_integrate_respa(int, int, int); void final_integrate_respa(int); double compute_scalar(); void write_restart(FILE *); diff --git a/src/fix_nve.cpp b/src/fix_nve.cpp index 53ea6f9395..a471c7217e 100644 --- a/src/fix_nve.cpp +++ b/src/fix_nve.cpp @@ -58,7 +58,7 @@ void FixNVE::init() allow for both per-type and per-atom mass ------------------------------------------------------------------------- */ -void FixNVE::initial_integrate() +void FixNVE::initial_integrate(int vflag) { double dtfm; @@ -137,14 +137,14 @@ void FixNVE::final_integrate() /* ---------------------------------------------------------------------- */ -void FixNVE::initial_integrate_respa(int ilevel, int flag) +void FixNVE::initial_integrate_respa(int vflag, int ilevel, int flag) { if (flag) return; // only used by NPT,NPH dtv = step_respa[ilevel]; dtf = 0.5 * step_respa[ilevel] * force->ftm2v; - if (ilevel == 0) initial_integrate(); + if (ilevel == 0) initial_integrate(vflag); else final_integrate(); } diff --git a/src/fix_nve.h b/src/fix_nve.h index 8c0e219bf2..88acbac073 100644 --- a/src/fix_nve.h +++ b/src/fix_nve.h @@ -23,9 +23,9 @@ class FixNVE : public Fix { FixNVE(class LAMMPS *, int, char **); int setmask(); virtual void init(); - virtual void initial_integrate(); + virtual void initial_integrate(int); virtual void final_integrate(); - void initial_integrate_respa(int, int); + void initial_integrate_respa(int, int, int); void final_integrate_respa(int); void reset_dt(); diff --git a/src/fix_nve_limit.cpp b/src/fix_nve_limit.cpp index 6195c4f944..d43684c708 100644 --- a/src/fix_nve_limit.cpp +++ b/src/fix_nve_limit.cpp @@ -67,7 +67,7 @@ void FixNVELimit::init() allow for both per-type and per-atom mass ------------------------------------------------------------------------- */ -void FixNVELimit::initial_integrate() +void FixNVELimit::initial_integrate(int vflag) { double dtfm,vsq,scale; @@ -184,14 +184,14 @@ void FixNVELimit::final_integrate() /* ---------------------------------------------------------------------- */ -void FixNVELimit::initial_integrate_respa(int ilevel, int flag) +void FixNVELimit::initial_integrate_respa(int vflag, int ilevel, int flag) { if (flag) return; // only used by NPT,NPH dtv = step_respa[ilevel]; dtf = 0.5 * step_respa[ilevel] * force->ftm2v; - if (ilevel == 0) initial_integrate(); + if (ilevel == 0) initial_integrate(vflag); else final_integrate(); } diff --git a/src/fix_nve_limit.h b/src/fix_nve_limit.h index a0c804ebf8..312300a76e 100644 --- a/src/fix_nve_limit.h +++ b/src/fix_nve_limit.h @@ -23,9 +23,9 @@ class FixNVELimit : public Fix { FixNVELimit(class LAMMPS *, int, char **); int setmask(); void init(); - void initial_integrate(); + void initial_integrate(int); void final_integrate(); - void initial_integrate_respa(int, int); + void initial_integrate_respa(int, int, int); void final_integrate_respa(int); void reset_dt(); double compute_scalar(); diff --git a/src/fix_nve_noforce.cpp b/src/fix_nve_noforce.cpp index 96da6dc9bf..2d8d91ec77 100644 --- a/src/fix_nve_noforce.cpp +++ b/src/fix_nve_noforce.cpp @@ -51,7 +51,7 @@ void FixNVENoforce::init() /* ---------------------------------------------------------------------- */ -void FixNVENoforce::initial_integrate() +void FixNVENoforce::initial_integrate(int vflag) { double **x = atom->x; double **v = atom->v; @@ -69,13 +69,13 @@ void FixNVENoforce::initial_integrate() /* ---------------------------------------------------------------------- */ -void FixNVENoforce::initial_integrate_respa(int ilevel, int flag) +void FixNVENoforce::initial_integrate_respa(int vflag, int ilevel, int flag) { if (flag) return; // only used by NPT,NPH dtv = step_respa[ilevel]; - if (ilevel == 0) initial_integrate(); + if (ilevel == 0) initial_integrate(vflag); } /* ---------------------------------------------------------------------- */ diff --git a/src/fix_nve_noforce.h b/src/fix_nve_noforce.h index 616a840373..3bbb0326be 100644 --- a/src/fix_nve_noforce.h +++ b/src/fix_nve_noforce.h @@ -23,8 +23,8 @@ class FixNVENoforce : public Fix { FixNVENoforce(class LAMMPS *, int, char **); int setmask(); void init(); - void initial_integrate(); - void initial_integrate_respa(int, int); + void initial_integrate(int); + void initial_integrate_respa(int, int, int); void reset_dt(); private: diff --git a/src/fix_nvt.cpp b/src/fix_nvt.cpp index cb194e3de1..f6807c127c 100644 --- a/src/fix_nvt.cpp +++ b/src/fix_nvt.cpp @@ -134,7 +134,7 @@ void FixNVT::init() /* ---------------------------------------------------------------------- */ -void FixNVT::setup() +void FixNVT::setup(int vflag) { t_target = t_start; // used by compute_scalar() t_current = temperature->compute_scalar(); @@ -142,7 +142,7 @@ void FixNVT::setup() /* ---------------------------------------------------------------------- */ -void FixNVT::initial_integrate() +void FixNVT::initial_integrate(int vflag) { double dtfm; @@ -218,7 +218,7 @@ void FixNVT::final_integrate() /* ---------------------------------------------------------------------- */ -void FixNVT::initial_integrate_respa(int ilevel, int flag) +void FixNVT::initial_integrate_respa(int vflag, int ilevel, int flag) { if (flag) return; // only used by NPT,NPH diff --git a/src/fix_nvt.h b/src/fix_nvt.h index 97618aa596..ff0ade8965 100644 --- a/src/fix_nvt.h +++ b/src/fix_nvt.h @@ -24,10 +24,10 @@ class FixNVT : public Fix { virtual ~FixNVT(); int setmask(); void init(); - void setup(); - virtual void initial_integrate(); + void setup(int); + virtual void initial_integrate(int); virtual void final_integrate(); - virtual void initial_integrate_respa(int,int); + virtual void initial_integrate_respa(int, int, int); void final_integrate_respa(int); double compute_scalar(); void write_restart(FILE *); diff --git a/src/fix_nvt_sllod.cpp b/src/fix_nvt_sllod.cpp index 2f9818622c..57db99d812 100644 --- a/src/fix_nvt_sllod.cpp +++ b/src/fix_nvt_sllod.cpp @@ -62,7 +62,7 @@ void FixNVTSlodd::init() /* ---------------------------------------------------------------------- */ -void FixNVTSlodd::initial_integrate() +void FixNVTSlodd::initial_integrate(int vflag) { double dtfm; @@ -191,7 +191,7 @@ void FixNVTSlodd::final_integrate() /* ---------------------------------------------------------------------- */ -void FixNVTSlodd::initial_integrate_respa(int ilevel, int flag) +void FixNVTSlodd::initial_integrate_respa(int vflag, int ilevel, int flag) { if (flag) return; // only used by NPT,NPH diff --git a/src/fix_nvt_sllod.h b/src/fix_nvt_sllod.h index abf96f8380..7adaac67eb 100644 --- a/src/fix_nvt_sllod.h +++ b/src/fix_nvt_sllod.h @@ -23,9 +23,9 @@ class FixNVTSlodd : public FixNVT { FixNVTSlodd(class LAMMPS *, int, char **); ~FixNVTSlodd() {} void init(); - void initial_integrate(); + void initial_integrate(int); void final_integrate(); - void initial_integrate_respa(int,int); + void initial_integrate_respa(int, int, int); }; } diff --git a/src/fix_orient_fcc.cpp b/src/fix_orient_fcc.cpp index 7599140b5c..b87edcc578 100644 --- a/src/fix_orient_fcc.cpp +++ b/src/fix_orient_fcc.cpp @@ -198,13 +198,13 @@ void FixOrientFCC::init_list(int id, NeighList *ptr) /* ---------------------------------------------------------------------- */ -void FixOrientFCC::setup() +void FixOrientFCC::setup(int vflag) { if (strcmp(update->integrate_style,"verlet") == 0) - post_force(1); + post_force(vflag); else { ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1); - post_force_respa(1,nlevels_respa-1,0); + post_force_respa(vflag,nlevels_respa-1,0); ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1); } } diff --git a/src/fix_orient_fcc.h b/src/fix_orient_fcc.h index 3293ba6700..2cf39634a9 100644 --- a/src/fix_orient_fcc.h +++ b/src/fix_orient_fcc.h @@ -42,7 +42,7 @@ class FixOrientFCC : public Fix { int setmask(); void init(); void init_list(int, class NeighList *); - void setup(); + void setup(int); void post_force(int); void post_force_respa(int, int, int); double compute_scalar(); diff --git a/src/fix_plane_force.cpp b/src/fix_plane_force.cpp index 104ff1162f..e4b7d98b72 100644 --- a/src/fix_plane_force.cpp +++ b/src/fix_plane_force.cpp @@ -45,15 +45,15 @@ int FixPlaneForce::setmask() /* ---------------------------------------------------------------------- */ -void FixPlaneForce::setup() +void FixPlaneForce::setup(int vflag) { if (strcmp(update->integrate_style,"verlet") == 0) - post_force(1); + post_force(vflag); else { int nlevels_respa = ((Respa *) update->integrate)->nlevels; for (int ilevel = 0; ilevel < nlevels_respa; ilevel++) { ((Respa *) update->integrate)->copy_flevel_f(ilevel); - post_force_respa(1,ilevel,0); + post_force_respa(vflag,ilevel,0); ((Respa *) update->integrate)->copy_f_flevel(ilevel); } } @@ -61,9 +61,9 @@ void FixPlaneForce::setup() /* ---------------------------------------------------------------------- */ -void FixPlaneForce::min_setup() +void FixPlaneForce::min_setup(int vflag) { - post_force(1); + post_force(vflag); } /* ---------------------------------------------------------------------- */ diff --git a/src/fix_plane_force.h b/src/fix_plane_force.h index 76c429d8b2..53cb5ba447 100644 --- a/src/fix_plane_force.h +++ b/src/fix_plane_force.h @@ -22,8 +22,8 @@ class FixPlaneForce : public Fix { public: FixPlaneForce(class LAMMPS *, int, char **); int setmask(); - void setup(); - void min_setup(); + void setup(int); + void min_setup(int); void post_force(int); void post_force_respa(int, int, int); void min_post_force(int); diff --git a/src/fix_rdf.cpp b/src/fix_rdf.cpp index ae981484da..1e94b53768 100644 --- a/src/fix_rdf.cpp +++ b/src/fix_rdf.cpp @@ -161,7 +161,7 @@ void FixRDF::init_list(int id, NeighList *ptr) /* ---------------------------------------------------------------------- */ -void FixRDF::setup() +void FixRDF::setup(int vflag) { if (first) end_of_step(); first = 0; diff --git a/src/fix_rdf.h b/src/fix_rdf.h index 4b2148402b..3ebbf54ec3 100644 --- a/src/fix_rdf.h +++ b/src/fix_rdf.h @@ -26,7 +26,7 @@ class FixRDF : public Fix { int setmask(); void init(); void init_list(int, class NeighList *); - void setup(); + void setup(int); void end_of_step(); private: diff --git a/src/fix_recenter.cpp b/src/fix_recenter.cpp index 604df0e436..e2f1406f86 100644 --- a/src/fix_recenter.cpp +++ b/src/fix_recenter.cpp @@ -132,7 +132,7 @@ void FixRecenter::init() /* ---------------------------------------------------------------------- */ -void FixRecenter::initial_integrate() +void FixRecenter::initial_integrate(int vflag) { // target COM // bounding box around domain works for both orthogonal and triclinic diff --git a/src/fix_recenter.h b/src/fix_recenter.h index b2a4e82656..38641bf90a 100644 --- a/src/fix_recenter.h +++ b/src/fix_recenter.h @@ -23,7 +23,7 @@ class FixRecenter : public Fix { FixRecenter(class LAMMPS *, int, char **); int setmask(); void init(); - void initial_integrate(); + void initial_integrate(int); private: int group2bit,scaleflag; diff --git a/src/fix_rigid.cpp b/src/fix_rigid.cpp index 4399d5043b..49ad4a0a72 100644 --- a/src/fix_rigid.cpp +++ b/src/fix_rigid.cpp @@ -219,10 +219,6 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) : if (screen) fprintf(screen,"%d rigid bodies with %d atoms\n",nbody,nsum); if (logfile) fprintf(logfile,"%d rigid bodies with %d atoms\n",nbody,nsum); } - - // zero fix_rigid virial in case pressure uses it before 1st fix_rigid call - - for (int n = 0; n < 6; n++) virial[n] = 0.0; } /* ---------------------------------------------------------------------- */ @@ -560,9 +556,9 @@ void FixRigid::init() /* ---------------------------------------------------------------------- */ -void FixRigid::setup() +void FixRigid::setup(int vflag) { - int i,ibody; + int i,n,ibody; // vcm = velocity of center-of-mass of each rigid body // fcm = force on center-of-mass of each rigid body @@ -660,21 +656,32 @@ void FixRigid::setup() torque[ibody][2] = all[ibody][5]; } + // virial setup before call to set_v + + if (vflag) v_setup(vflag); + else evflag = 0; + // set velocities from angmom & omega - // guestimate virial as 2x the set_v contribution for (ibody = 0; ibody < nbody; ibody++) omega_from_mq(angmom[ibody],ex_space[ibody],ey_space[ibody], ez_space[ibody],inertia[ibody],omega[ibody]); + set_v(); - for (int n = 0; n < 6; n++) virial[n] = 0.0; - set_v(1); - for (int n = 0; n < 6; n++) virial[n] *= 2.0; + // guestimate virial as 2x the set_v contribution + + if (vflag_global) + for (n = 0; n < 6; n++) virial[n] *= 2.0; + if (vflag_atom) { + for (i = 0; i < nlocal; i++) + for (n = 0; n < 6; n++) + vatom[i][n] *= 2.0; + } } /* ---------------------------------------------------------------------- */ -void FixRigid::initial_integrate() +void FixRigid::initial_integrate(int vflag) { double dtfm; @@ -708,12 +715,15 @@ void FixRigid::initial_integrate() ex_space[ibody],ey_space[ibody],ez_space[ibody]); } + // virial setup before call to set_xv + + if (vflag) v_setup(vflag); + else evflag = 0; + // set coords and velocities if atoms in rigid bodies // from quarternion and omega - int vflag = 0; - if (pressure_flag || output->next_thermo == update->ntimestep) vflag = 1; - set_xv(vflag); + set_xv(); } /* ---------------------------------------------------------------------- */ @@ -862,19 +872,18 @@ void FixRigid::final_integrate() } // set velocities from angmom & omega - + // virial is already setup from initial_integrate + for (ibody = 0; ibody < nbody; ibody++) omega_from_mq(angmom[ibody],ex_space[ibody],ey_space[ibody], ez_space[ibody],inertia[ibody],omega[ibody]); - int vflag = 0; - if (pressure_flag || output->next_thermo == update->ntimestep) vflag = 1; - set_v(vflag); + set_v(); } /* ---------------------------------------------------------------------- */ -void FixRigid::initial_integrate_respa(int ilevel, int flag) +void FixRigid::initial_integrate_respa(int vflag, int ilevel, int flag) { if (flag) return; // only used by NPT,NPH @@ -882,7 +891,7 @@ void FixRigid::initial_integrate_respa(int ilevel, int flag) dtf = 0.5 * step_respa[ilevel] * force->ftm2v; dtq = 0.5 * step_respa[ilevel]; - if (ilevel == 0) initial_integrate(); + if (ilevel == 0) initial_integrate(vflag); else final_integrate(); } @@ -1220,12 +1229,13 @@ void FixRigid::exyz_from_q(double *q, double *ex, double *ey, double *ez) v = Vcm + (W cross (x - Xcm)) ------------------------------------------------------------------------- */ -void FixRigid::set_xv(int vflag) +void FixRigid::set_xv() { int ibody; int xbox,ybox,zbox; - double vold0,vold1,vold2,fc0,fc1,fc2,massone,x0,x1,x2; + double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone; double xy,xz,yz; + double vr[6]; int *image = atom->image; double **x = atom->x; @@ -1246,10 +1256,6 @@ void FixRigid::set_xv(int vflag) double **f = atom->f; int *type = atom->type; - // zero out fix_rigid virial - - if (vflag) for (int n = 0; n < 6; n++) virial[n] = 0.0; - for (int i = 0; i < nlocal; i++) { if (body[i] < 0) continue; ibody = body[i]; @@ -1258,16 +1264,21 @@ void FixRigid::set_xv(int vflag) ybox = (image[i] >> 10 & 1023) - 512; zbox = (image[i] >> 20) - 512; - // save old positions and velocities for virial contribution + // save old positions and velocities for virial - if (vflag) { - x0 = x[i][0] + xbox*xprd; - x1 = x[i][1] + ybox*yprd; - x2 = x[i][2] + zbox*zprd; - - vold0 = v[i][0]; - vold1 = v[i][1]; - vold2 = v[i][2]; + if (evflag) { + if (triclinic == 0) { + x0 = x[i][0] + xbox*xprd; + x1 = x[i][1] + ybox*yprd; + x2 = x[i][2] + zbox*zprd; + } else { + x0 = x[i][0] + xbox*xprd + ybox*xy + zbox*xz; + x1 = x[i][1] + ybox*yprd + zbox*yz; + x2 = x[i][2] + zbox*zprd; + } + v0 = v[i][0]; + v1 = v[i][1]; + v2 = v[i][2]; } // x = displacement from center-of-mass, based on body orientation @@ -1304,20 +1315,26 @@ void FixRigid::set_xv(int vflag) x[i][2] += xcm[ibody][2] - zbox*zprd; } - // compute body constraint forces for virial + // virial = unwrapped coords dotted into body constraint force + // body constraint force = implied force due to v change minus f external + // assume f does not include forces internal to body + // 1/2 factor b/c final_integrate contributes other half + // assume per-atom contribution is due to constraint force on that atom - if (vflag) { + if (evflag) { massone = mass[type[i]]; - fc0 = massone*(v[i][0] - vold0)/dtf - f[i][0]; - fc1 = massone*(v[i][1] - vold1)/dtf - f[i][1]; - fc2 = massone*(v[i][2] - vold2)/dtf - f[i][2]; + fc0 = massone*(v[i][0] - v0)/dtf - f[i][0]; + fc1 = massone*(v[i][1] - v1)/dtf - f[i][1]; + fc2 = massone*(v[i][2] - v2)/dtf - f[i][2]; - virial[0] += 0.5*fc0*x0; - virial[1] += 0.5*fc1*x1; - virial[2] += 0.5*fc2*x2; - virial[3] += 0.5*fc1*x0; - virial[4] += 0.5*fc2*x0; - virial[5] += 0.5*fc2*x1; + vr[0] = 0.5*x0*fc0; + vr[1] = 0.5*x1*fc1; + vr[2] = 0.5*x2*fc2; + vr[3] = 0.5*x0*fc1; + vr[4] = 0.5*x0*fc2; + vr[5] = 0.5*x1*fc2; + + v_tally(1,&i,1.0,vr); } } } @@ -1327,22 +1344,22 @@ void FixRigid::set_xv(int vflag) v = Vcm + (W cross (x - Xcm)) ------------------------------------------------------------------------- */ -void FixRigid::set_v(int vflag) +void FixRigid::set_v() { - double **v = atom->v; - int nlocal = atom->nlocal; - int ibody; + int xbox,ybox,zbox; double xunwrap,yunwrap,zunwrap,dx,dy,dz; - double vold0,vold1,vold2,fc0,fc1,fc2,massone; + double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone; double xy,xz,yz; + double vr[6]; double *mass = atom->mass; double **f = atom->f; + double **v = atom->v; double **x = atom->x; int *type = atom->type; int *image = atom->image; - int xbox,ybox,zbox; + int nlocal = atom->nlocal; double xprd = domain->xprd; double yprd = domain->yprd; @@ -1369,45 +1386,50 @@ void FixRigid::set_v(int vflag) // save old velocities for virial - if (vflag) { - vold0 = v[i][0]; - vold1 = v[i][1]; - vold2 = v[i][2]; + if (evflag) { + v0 = v[i][0]; + v1 = v[i][1]; + v2 = v[i][2]; } v[i][0] = omega[ibody][1]*dz - omega[ibody][2]*dy + vcm[ibody][0]; v[i][1] = omega[ibody][2]*dx - omega[ibody][0]*dz + vcm[ibody][1]; v[i][2] = omega[ibody][0]*dy - omega[ibody][1]*dx + vcm[ibody][2]; - // compute body constraint forces for virial - // use unwrapped atom positions + // virial = unwrapped coords dotted into body constraint force + // body constraint force = implied force due to v change minus f external + // assume f does not include forces internal to body + // 1/2 factor b/c initial_integrate contributes other half + // assume per-atom contribution is due to constraint force on that atom - if (vflag) { + if (evflag) { massone = mass[type[i]]; - fc0 = massone*(v[i][0] - vold0)/dtf - f[i][0]; - fc1 = massone*(v[i][1] - vold1)/dtf - f[i][1]; - fc2 = massone*(v[i][2] - vold2)/dtf - f[i][2]; + fc0 = massone*(v[i][0] - v0)/dtf - f[i][0]; + fc1 = massone*(v[i][1] - v1)/dtf - f[i][1]; + fc2 = massone*(v[i][2] - v2)/dtf - f[i][2]; xbox = (image[i] & 1023) - 512; ybox = (image[i] >> 10 & 1023) - 512; zbox = (image[i] >> 20) - 512; if (triclinic == 0) { - xunwrap = x[i][0] + xbox*xprd; - yunwrap = x[i][1] + ybox*yprd; - zunwrap = x[i][2] + zbox*zprd; + x0 = x[i][0] + xbox*xprd; + x1 = x[i][1] + ybox*yprd; + x2 = x[i][2] + zbox*zprd; } else { - xunwrap = x[i][0] + xbox*xprd + ybox*xy + zbox*xz; - yunwrap = x[i][1] + ybox*yprd + zbox*yz; - zunwrap = x[i][2] + zbox*zprd; + x0 = x[i][0] + xbox*xprd + ybox*xy + zbox*xz; + x1 = x[i][1] + ybox*yprd + zbox*yz; + x2 = x[i][2] + zbox*zprd; } - virial[0] += 0.5*fc0*xunwrap; - virial[1] += 0.5*fc1*yunwrap; - virial[2] += 0.5*fc2*zunwrap; - virial[3] += 0.5*fc1*xunwrap; - virial[4] += 0.5*fc2*xunwrap; - virial[5] += 0.5*fc2*yunwrap; + vr[0] = 0.5*x0*fc0; + vr[1] = 0.5*x1*fc1; + vr[2] = 0.5*x2*fc2; + vr[3] = 0.5*x0*fc1; + vr[4] = 0.5*x0*fc2; + vr[5] = 0.5*x1*fc2; + + v_tally(1,&i,1.0,vr); } } } @@ -1421,6 +1443,7 @@ double FixRigid::memory_usage() int nmax = atom->nmax; double bytes = nmax * sizeof(int); bytes += nmax*3 * sizeof(double); + bytes += maxvatom*6 * sizeof(double); return bytes; } diff --git a/src/fix_rigid.h b/src/fix_rigid.h index 4b59512022..46d26a4fb5 100644 --- a/src/fix_rigid.h +++ b/src/fix_rigid.h @@ -24,10 +24,10 @@ class FixRigid : public Fix { ~FixRigid(); int setmask(); void init(); - void setup(); - void initial_integrate(); + void setup(int); + void initial_integrate(int); void final_integrate(); - void initial_integrate_respa(int, int); + void initial_integrate_respa(int, int, int); void final_integrate_respa(int); double memory_usage(); @@ -77,8 +77,8 @@ class FixRigid : public Fix { void omega_from_mq(double *, double *, double *, double *, double *, double *); void exyz_from_q(double *, double *, double *, double *); - void set_xv(int); - void set_v(int); + void set_xv(); + void set_v(); }; } diff --git a/src/fix_set_force.cpp b/src/fix_set_force.cpp index 795723959c..70c0666634 100644 --- a/src/fix_set_force.cpp +++ b/src/fix_set_force.cpp @@ -66,23 +66,23 @@ void FixSetForce::init() /* ---------------------------------------------------------------------- */ -void FixSetForce::setup() +void FixSetForce::setup(int vflag) { if (strcmp(update->integrate_style,"verlet") == 0) - post_force(1); + post_force(vflag); else for (int ilevel = 0; ilevel < nlevels_respa; ilevel++) { ((Respa *) update->integrate)->copy_flevel_f(ilevel); - post_force_respa(1,ilevel,0); + post_force_respa(vflag,ilevel,0); ((Respa *) update->integrate)->copy_f_flevel(ilevel); } } /* ---------------------------------------------------------------------- */ -void FixSetForce::min_setup() +void FixSetForce::min_setup(int vflag) { - post_force(1); + post_force(vflag); } /* ---------------------------------------------------------------------- */ diff --git a/src/fix_set_force.h b/src/fix_set_force.h index 15ea1a3b70..5f0f6d88b7 100644 --- a/src/fix_set_force.h +++ b/src/fix_set_force.h @@ -23,8 +23,8 @@ class FixSetForce : public Fix { FixSetForce(class LAMMPS *, int, char **); int setmask(); void init(); - void setup(); - void min_setup(); + void setup(int); + void min_setup(int); void post_force(int); void post_force_respa(int, int, int); void min_post_force(int); diff --git a/src/fix_shake.cpp b/src/fix_shake.cpp index d598dc8db8..ecc2c78609 100644 --- a/src/fix_shake.cpp +++ b/src/fix_shake.cpp @@ -379,7 +379,7 @@ void FixShake::init() SHAKE as pre-integrator constraint ------------------------------------------------------------------------- */ -void FixShake::setup() +void FixShake::setup(int vflag) { pre_neighbor(); @@ -398,14 +398,14 @@ void FixShake::setup() if (strcmp(update->integrate_style,"verlet") == 0) { dtv = update->dt; dtfsq = 0.5 * update->dt * update->dt * force->ftm2v; - post_force(1); + post_force(vflag); dtfsq = update->dt * update->dt * force->ftm2v; } else { 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(1,nlevels_respa-1,0); + post_force_respa(vflag,nlevels_respa-1,0); ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1); dtf_inner = step_respa[0] * force->ftm2v; } @@ -490,7 +490,7 @@ void FixShake::pre_neighbor() compute the force adjustment for SHAKE constraint ------------------------------------------------------------------------- */ -void FixShake::post_force(int vflag_in) +void FixShake::post_force(int vflag) { if (update->ntimestep == next_output) stats(); @@ -502,10 +502,10 @@ void FixShake::post_force(int vflag_in) if (nprocs > 1) comm->comm_fix(this); - // zero out SHAKE contribution to virial + // virial setup - vflag = vflag_in; - if (vflag) for (int n = 0; n < 6; n++) virial[n] = 0.0; + if (vflag) v_setup(vflag); + else evflag = 0; // loop over clusters @@ -1148,6 +1148,9 @@ void FixShake::unconstrained_update() void FixShake::shake2(int m) { + int nlist,list[2]; + double v[6]; + // local atom IDs and constraint distances int i0 = atom->map(shake_atom[m][0]); @@ -1218,18 +1221,19 @@ void FixShake::shake2(int m) f[i1][2] -= lamda*r01[2]; } - if (vflag) { - int factor = 0; - if (i0 < nlocal) factor++; - if (i1 < nlocal) factor++; - double rfactor = 0.5 * factor; + if (evflag) { + nlist = 0; + if (i0 < nlocal) list[nlist++] = i0; + if (i1 < nlocal) list[nlist++] = i1; - virial[0] += rfactor*lamda*r01[0]*r01[0]; - virial[1] += rfactor*lamda*r01[1]*r01[1]; - virial[2] += rfactor*lamda*r01[2]*r01[2]; - virial[3] += rfactor*lamda*r01[0]*r01[1]; - virial[4] += rfactor*lamda*r01[0]*r01[2]; - virial[5] += rfactor*lamda*r01[1]*r01[2]; + v[0] = lamda*r01[0]*r01[0]; + v[1] = lamda*r01[1]*r01[1]; + v[2] = lamda*r01[2]*r01[2]; + v[3] = lamda*r01[0]*r01[1]; + v[4] = lamda*r01[0]*r01[2]; + v[5] = lamda*r01[1]*r01[2]; + + v_tally(nlist,list,2.0,v); } } @@ -1237,6 +1241,9 @@ void FixShake::shake2(int m) void FixShake::shake3(int m) { + int nlist,list[3]; + double v[6]; + // local atom IDs and constraint distances int i0 = atom->map(shake_atom[m][0]); @@ -1371,26 +1378,20 @@ void FixShake::shake3(int m) f[i2][2] -= lamda02*r02[2]; } - if (vflag) { - int factor = 0; - if (i0 < nlocal) factor++; - if (i1 < nlocal) factor++; - if (i2 < nlocal) factor++; - double rfactor = factor/3.0; + if (evflag) { + nlist = 0; + if (i0 < nlocal) list[nlist++] = i0; + if (i1 < nlocal) list[nlist++] = i1; + if (i2 < nlocal) list[nlist++] = i2; - virial[0] += rfactor*lamda01*r01[0]*r01[0]; - virial[1] += rfactor*lamda01*r01[1]*r01[1]; - virial[2] += rfactor*lamda01*r01[2]*r01[2]; - virial[3] += rfactor*lamda01*r01[0]*r01[1]; - virial[4] += rfactor*lamda01*r01[0]*r01[2]; - virial[5] += rfactor*lamda01*r01[1]*r01[2]; + v[0] = lamda01*r01[0]*r01[0] + lamda02*r02[0]*r02[0]; + v[1] = lamda01*r01[1]*r01[1] + lamda02*r02[1]*r02[1]; + v[2] = lamda01*r01[2]*r01[2] + lamda02*r02[2]*r02[2]; + v[3] = lamda01*r01[0]*r01[1] + lamda02*r02[0]*r02[1]; + v[4] = lamda01*r01[0]*r01[2] + lamda02*r02[0]*r02[2]; + v[5] = lamda01*r01[1]*r01[2] + lamda02*r02[1]*r02[2]; - virial[0] += rfactor*lamda02*r02[0]*r02[0]; - virial[1] += rfactor*lamda02*r02[1]*r02[1]; - virial[2] += rfactor*lamda02*r02[2]*r02[2]; - virial[3] += rfactor*lamda02*r02[0]*r02[1]; - virial[4] += rfactor*lamda02*r02[0]*r02[2]; - virial[5] += rfactor*lamda02*r02[1]*r02[2]; + v_tally(nlist,list,3.0,v); } } @@ -1398,6 +1399,9 @@ void FixShake::shake3(int m) void FixShake::shake4(int m) { + int nlist,list[4]; + double v[6]; + // local atom IDs and constraint distances int i0 = atom->map(shake_atom[m][0]); @@ -1608,34 +1612,21 @@ void FixShake::shake4(int m) f[i3][2] -= lamda03*r03[2]; } - if (vflag) { - int factor = 0; - if (i0 < nlocal) factor++; - if (i1 < nlocal) factor++; - if (i2 < nlocal) factor++; - if (i3 < nlocal) factor++; - double rfactor = 0.25*factor; + if (evflag) { + nlist = 0; + if (i0 < nlocal) list[nlist++] = i0; + if (i1 < nlocal) list[nlist++] = i1; + if (i2 < nlocal) list[nlist++] = i2; + if (i3 < nlocal) list[nlist++] = i3; - virial[0] += rfactor*lamda01*r01[0]*r01[0]; - virial[1] += rfactor*lamda01*r01[1]*r01[1]; - virial[2] += rfactor*lamda01*r01[2]*r01[2]; - virial[3] += rfactor*lamda01*r01[0]*r01[1]; - virial[4] += rfactor*lamda01*r01[0]*r01[2]; - virial[5] += rfactor*lamda01*r01[1]*r01[2]; + v[0] = lamda01*r01[0]*r01[0]+lamda02*r02[0]*r02[0]+lamda03*r03[0]*r03[0]; + v[1] = lamda01*r01[1]*r01[1]+lamda02*r02[1]*r02[1]+lamda03*r03[1]*r03[1]; + v[2] = lamda01*r01[2]*r01[2]+lamda02*r02[2]*r02[2]+lamda03*r03[2]*r03[2]; + v[3] = lamda01*r01[0]*r01[1]+lamda02*r02[0]*r02[1]+lamda03*r03[0]*r03[1]; + v[4] = lamda01*r01[0]*r01[2]+lamda02*r02[0]*r02[2]+lamda03*r03[0]*r03[2]; + v[5] = lamda01*r01[1]*r01[2]+lamda02*r02[1]*r02[2]+lamda03*r03[1]*r03[2]; - virial[0] += rfactor*lamda02*r02[0]*r02[0]; - virial[1] += rfactor*lamda02*r02[1]*r02[1]; - virial[2] += rfactor*lamda02*r02[2]*r02[2]; - virial[3] += rfactor*lamda02*r02[0]*r02[1]; - virial[4] += rfactor*lamda02*r02[0]*r02[2]; - virial[5] += rfactor*lamda02*r02[1]*r02[2]; - - virial[0] += rfactor*lamda03*r03[0]*r03[0]; - virial[1] += rfactor*lamda03*r03[1]*r03[1]; - virial[2] += rfactor*lamda03*r03[2]*r03[2]; - virial[3] += rfactor*lamda03*r03[0]*r03[1]; - virial[4] += rfactor*lamda03*r03[0]*r03[2]; - virial[5] += rfactor*lamda03*r03[1]*r03[2]; + v_tally(nlist,list,4.0,v); } } @@ -1643,6 +1634,9 @@ void FixShake::shake4(int m) void FixShake::shake3angle(int m) { + int nlist,list[3]; + double v[6]; + // local atom IDs and constraint distances int i0 = atom->map(shake_atom[m][0]); @@ -1845,33 +1839,20 @@ void FixShake::shake3angle(int m) f[i2][2] -= lamda02*r02[2] + lamda12*r12[2]; } - if (vflag) { - int factor = 0; - if (i0 < nlocal) factor++; - if (i1 < nlocal) factor++; - if (i2 < nlocal) factor++; - double rfactor = factor/3.0; + if (evflag) { + nlist = 0; + if (i0 < nlocal) list[nlist++] = i0; + if (i1 < nlocal) list[nlist++] = i1; + if (i2 < nlocal) list[nlist++] = i2; - virial[0] += rfactor*lamda01*r01[0]*r01[0]; - virial[1] += rfactor*lamda01*r01[1]*r01[1]; - virial[2] += rfactor*lamda01*r01[2]*r01[2]; - virial[3] += rfactor*lamda01*r01[0]*r01[1]; - virial[4] += rfactor*lamda01*r01[0]*r01[2]; - virial[5] += rfactor*lamda01*r01[1]*r01[2]; + v[0] = lamda01*r01[0]*r01[0]+lamda02*r02[0]*r02[0]+lamda12*r12[0]*r12[0]; + v[1] = lamda01*r01[1]*r01[1]+lamda02*r02[1]*r02[1]+lamda12*r12[1]*r12[1]; + v[2] = lamda01*r01[2]*r01[2]+lamda02*r02[2]*r02[2]+lamda12*r12[2]*r12[2]; + v[3] = lamda01*r01[0]*r01[1]+lamda02*r02[0]*r02[1]+lamda12*r12[0]*r12[1]; + v[4] = lamda01*r01[0]*r01[2]+lamda02*r02[0]*r02[2]+lamda12*r12[0]*r12[2]; + v[5] = lamda01*r01[1]*r01[2]+lamda02*r02[1]*r02[2]+lamda12*r12[1]*r12[2]; - virial[0] += rfactor*lamda02*r02[0]*r02[0]; - virial[1] += rfactor*lamda02*r02[1]*r02[1]; - virial[2] += rfactor*lamda02*r02[2]*r02[2]; - virial[3] += rfactor*lamda02*r02[0]*r02[1]; - virial[4] += rfactor*lamda02*r02[0]*r02[2]; - virial[5] += rfactor*lamda02*r02[1]*r02[2]; - - virial[0] += rfactor*lamda12*r12[0]*r12[0]; - virial[1] += rfactor*lamda12*r12[1]*r12[1]; - virial[2] += rfactor*lamda12*r12[2]*r12[2]; - virial[3] += rfactor*lamda12*r12[0]*r12[1]; - virial[4] += rfactor*lamda12*r12[0]*r12[2]; - virial[5] += rfactor*lamda12*r12[1]*r12[2]; + v_tally(nlist,list,3.0,v); } } @@ -2064,6 +2045,7 @@ double FixShake::memory_usage() bytes += nmax*4 * sizeof(int); bytes += nmax*3 * sizeof(int); bytes += nmax*3 * sizeof(double); + bytes += maxvatom*6 * sizeof(double); return bytes; } @@ -2199,7 +2181,7 @@ int FixShake::unpack_exchange(int nlocal, double *buf) rRESPA updating of atom coords is done with full v, but only portions of f ------------------------------------------------------------------------- */ -void FixShake::post_force_respa(int vflag_in, int ilevel, int iloop) +void FixShake::post_force_respa(int vflag, int ilevel, int iloop) { // call stats only on outermost level @@ -2242,10 +2224,10 @@ void FixShake::post_force_respa(int vflag_in, int ilevel, int iloop) if (nprocs > 1) comm->comm_fix(this); - // zero out SHAKE contribution to virial + // virial setup - vflag = vflag_in; - if (vflag) for (int n = 0; n < 6; n++) virial[n] = 0.0; + if (vflag) v_setup(vflag); + else evflag = 0; // loop over clusters diff --git a/src/fix_shake.h b/src/fix_shake.h index 73c2ed4689..3512bc43ea 100644 --- a/src/fix_shake.h +++ b/src/fix_shake.h @@ -24,7 +24,7 @@ class FixShake : public Fix { ~FixShake(); int setmask(); void init(); - void setup(); + void setup(int); void pre_neighbor(); void post_force(int); void post_force_respa(int, int, int); diff --git a/src/fix_spring.cpp b/src/fix_spring.cpp index 771f8f81c3..479a639327 100644 --- a/src/fix_spring.cpp +++ b/src/fix_spring.cpp @@ -106,13 +106,13 @@ void FixSpring::init() /* ---------------------------------------------------------------------- */ -void FixSpring::setup() +void FixSpring::setup(int vflag) { if (strcmp(update->integrate_style,"verlet") == 0) - post_force(1); + post_force(vflag); else { ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1); - post_force_respa(1,nlevels_respa-1,0); + post_force_respa(vflag,nlevels_respa-1,0); ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1); } } diff --git a/src/fix_spring.h b/src/fix_spring.h index 66e917cf23..381e221d26 100644 --- a/src/fix_spring.h +++ b/src/fix_spring.h @@ -23,7 +23,7 @@ class FixSpring : public Fix { FixSpring(class LAMMPS *, int, char **); int setmask(); void init(); - void setup(); + void setup(int); void post_force(int); void post_force_respa(int, int, int); double compute_vector(int); diff --git a/src/fix_spring_rg.cpp b/src/fix_spring_rg.cpp index 0588297b7e..6337967ddb 100644 --- a/src/fix_spring_rg.cpp +++ b/src/fix_spring_rg.cpp @@ -74,13 +74,13 @@ void FixSpringRG::init() /* ---------------------------------------------------------------------- */ -void FixSpringRG::setup() +void FixSpringRG::setup(int vflag) { if (strcmp(update->integrate_style,"verlet") == 0) - post_force(1); + post_force(vflag); else { ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1); - post_force_respa(1,nlevels_respa-1,0); + post_force_respa(vflag,nlevels_respa-1,0); ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1); } } diff --git a/src/fix_spring_rg.h b/src/fix_spring_rg.h index ab91c57aa7..a922a41500 100644 --- a/src/fix_spring_rg.h +++ b/src/fix_spring_rg.h @@ -23,7 +23,7 @@ class FixSpringRG : public Fix { FixSpringRG(class LAMMPS *, int, char **); int setmask(); void init(); - void setup(); + void setup(int); void post_force(int); void post_force_respa(int, int, int); diff --git a/src/fix_spring_self.cpp b/src/fix_spring_self.cpp index f498c5935d..1e84ca31b4 100644 --- a/src/fix_spring_self.cpp +++ b/src/fix_spring_self.cpp @@ -105,13 +105,13 @@ void FixSpringSelf::init() /* ---------------------------------------------------------------------- */ -void FixSpringSelf::setup() +void FixSpringSelf::setup(int vflag) { if (strcmp(update->integrate_style,"verlet") == 0) - post_force(1); + post_force(vflag); else { ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1); - post_force_respa(1,nlevels_respa-1,0); + post_force_respa(vflag,nlevels_respa-1,0); ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1); } } diff --git a/src/fix_spring_self.h b/src/fix_spring_self.h index e098ff3b87..7eed32e10b 100644 --- a/src/fix_spring_self.h +++ b/src/fix_spring_self.h @@ -24,7 +24,7 @@ class FixSpringSelf : public Fix { ~FixSpringSelf(); int setmask(); void init(); - void setup(); + void setup(int); void post_force(int); void post_force_respa(int, int, int); diff --git a/src/fix_tmd.cpp b/src/fix_tmd.cpp index 96213bc2e2..5b41637ba8 100644 --- a/src/fix_tmd.cpp +++ b/src/fix_tmd.cpp @@ -180,7 +180,7 @@ void FixTMD::init() /* ---------------------------------------------------------------------- */ -void FixTMD::initial_integrate() +void FixTMD::initial_integrate(int vflag) { double a,b,c,d,e; double dx,dy,dz,dxkt,dykt,dzkt; @@ -321,14 +321,14 @@ void FixTMD::initial_integrate() /* ---------------------------------------------------------------------- */ -void FixTMD::initial_integrate_respa(int ilevel, int flag) +void FixTMD::initial_integrate_respa(int vflag, int ilevel, int flag) { if (flag) return; // only used by NPT,NPH dtv = step_respa[ilevel]; dtf = step_respa[ilevel] * force->ftm2v; - if (ilevel == 0) initial_integrate(); + if (ilevel == 0) initial_integrate(vflag); } /* ---------------------------------------------------------------------- diff --git a/src/fix_tmd.h b/src/fix_tmd.h index 78f61633ca..7ea5582fed 100644 --- a/src/fix_tmd.h +++ b/src/fix_tmd.h @@ -25,8 +25,8 @@ class FixTMD : public Fix { ~FixTMD(); int setmask(); void init(); - void initial_integrate(); - void initial_integrate_respa(int,int); + void initial_integrate(int); + void initial_integrate_respa(int, int, int); double memory_usage(); void grow_arrays(int); diff --git a/src/fix_viscous.cpp b/src/fix_viscous.cpp index 8fdbb0da89..b6c8b18db8 100644 --- a/src/fix_viscous.cpp +++ b/src/fix_viscous.cpp @@ -76,13 +76,13 @@ void FixViscous::init() /* ---------------------------------------------------------------------- */ -void FixViscous::setup() +void FixViscous::setup(int vflag) { if (strcmp(update->integrate_style,"verlet") == 0) - post_force(1); + post_force(vflag); else { ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1); - post_force_respa(1,nlevels_respa-1,0); + post_force_respa(vflag,nlevels_respa-1,0); ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1); } } diff --git a/src/fix_viscous.h b/src/fix_viscous.h index a682c99bde..eddf49d735 100644 --- a/src/fix_viscous.h +++ b/src/fix_viscous.h @@ -24,7 +24,7 @@ class FixViscous : public Fix { ~FixViscous(); int setmask(); void init(); - void setup(); + void setup(int); void post_force(int); void post_force_respa(int, int, int); diff --git a/src/fix_wall_lj126.cpp b/src/fix_wall_lj126.cpp index e7fa0af8fc..0755602aac 100644 --- a/src/fix_wall_lj126.cpp +++ b/src/fix_wall_lj126.cpp @@ -109,22 +109,22 @@ void FixWallLJ126::init() /* ---------------------------------------------------------------------- */ -void FixWallLJ126::setup() +void FixWallLJ126::setup(int vflag) { if (strcmp(update->integrate_style,"verlet") == 0) - post_force(1); + post_force(vflag); else { ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1); - post_force_respa(1,nlevels_respa-1,0); + post_force_respa(vflag,nlevels_respa-1,0); ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1); } } /* ---------------------------------------------------------------------- */ -void FixWallLJ126::min_setup() +void FixWallLJ126::min_setup(int vflag) { - post_force(1); + post_force(vflag); } /* ---------------------------------------------------------------------- */ diff --git a/src/fix_wall_lj126.h b/src/fix_wall_lj126.h index 2118af4f43..39678cc8fd 100644 --- a/src/fix_wall_lj126.h +++ b/src/fix_wall_lj126.h @@ -23,8 +23,8 @@ class FixWallLJ126 : public Fix { FixWallLJ126(class LAMMPS *, int, char **); int setmask(); void init(); - void setup(); - void min_setup(); + void setup(int); + void min_setup(int); void post_force(int); void post_force_respa(int, int, int); void min_post_force(int); diff --git a/src/fix_wall_lj93.cpp b/src/fix_wall_lj93.cpp index 4c433d7f91..34cb4fd09b 100644 --- a/src/fix_wall_lj93.cpp +++ b/src/fix_wall_lj93.cpp @@ -106,22 +106,22 @@ void FixWallLJ93::init() /* ---------------------------------------------------------------------- */ -void FixWallLJ93::setup() +void FixWallLJ93::setup(int vflag) { if (strcmp(update->integrate_style,"verlet") == 0) - post_force(1); + post_force(vflag); else { ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1); - post_force_respa(1,nlevels_respa-1,0); + post_force_respa(vflag,nlevels_respa-1,0); ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1); } } /* ---------------------------------------------------------------------- */ -void FixWallLJ93::min_setup() +void FixWallLJ93::min_setup(int vflag) { - post_force(1); + post_force(vflag); } /* ---------------------------------------------------------------------- */ diff --git a/src/fix_wall_lj93.h b/src/fix_wall_lj93.h index 9a6d2cc6f5..b5a5986086 100644 --- a/src/fix_wall_lj93.h +++ b/src/fix_wall_lj93.h @@ -24,8 +24,8 @@ class FixWallLJ93 : public Fix { ~FixWallLJ93() {} int setmask(); void init(); - void setup(); - void min_setup(); + void setup(int); + void min_setup(int); void post_force(int); void post_force_respa(int, int, int); void min_post_force(int); diff --git a/src/fix_wall_reflect.cpp b/src/fix_wall_reflect.cpp index dbb4e260a5..a1aa470f9b 100644 --- a/src/fix_wall_reflect.cpp +++ b/src/fix_wall_reflect.cpp @@ -74,7 +74,7 @@ void FixWallReflect::init() /* ---------------------------------------------------------------------- */ -void FixWallReflect::initial_integrate() +void FixWallReflect::initial_integrate(int vflag) { double xlo = domain->boxlo[0]; double xhi = domain->boxhi[0]; diff --git a/src/fix_wall_reflect.h b/src/fix_wall_reflect.h index fecf952ab4..53c588ad72 100644 --- a/src/fix_wall_reflect.h +++ b/src/fix_wall_reflect.h @@ -23,7 +23,7 @@ class FixWallReflect : public Fix { FixWallReflect(class LAMMPS *, int, char **); int setmask(); void init(); - void initial_integrate(); + void initial_integrate(int); private: int xloflag,xhiflag,yloflag,yhiflag,zloflag,zhiflag; diff --git a/src/min_cg.cpp b/src/min_cg.cpp index 72bafad0b4..853791ef39 100644 --- a/src/min_cg.cpp +++ b/src/min_cg.cpp @@ -259,7 +259,7 @@ void MinCG::setup() if (force->newton) comm->reverse_communicate(); - modify->setup(); + modify->setup(vflag); output->setup(1); } diff --git a/src/modify.cpp b/src/modify.cpp index 2c2dea521e..20ab5e9324 100644 --- a/src/modify.cpp +++ b/src/modify.cpp @@ -171,22 +171,22 @@ void Modify::init() setup for run, calls setup() of all fixes ------------------------------------------------------------------------- */ -void Modify::setup() +void Modify::setup(int vflag) { if (update->whichflag == 0) - for (int i = 0; i < nfix; i++) fix[i]->setup(); + for (int i = 0; i < nfix; i++) fix[i]->setup(vflag); else - for (int i = 0; i < nfix; i++) fix[i]->min_setup(); + for (int i = 0; i < nfix; i++) fix[i]->min_setup(vflag); } /* ---------------------------------------------------------------------- 1st half of integrate call only for relevant fixes ------------------------------------------------------------------------- */ -void Modify::initial_integrate() +void Modify::initial_integrate(int vflag) { for (int i = 0; i < n_initial_integrate; i++) - fix[list_initial_integrate[i]]->initial_integrate(); + fix[list_initial_integrate[i]]->initial_integrate(vflag); } /* ---------------------------------------------------------------------- @@ -259,10 +259,11 @@ double Modify::thermo_energy() 1st half of rRESPA integrate call only for relevant fixes ------------------------------------------------------------------------- */ -void Modify::initial_integrate_respa(int ilevel, int flag) +void Modify::initial_integrate_respa(int vflag, int ilevel, int flag) { for (int i = 0; i < n_initial_integrate_respa; i++) - fix[list_initial_integrate_respa[i]]->initial_integrate_respa(ilevel,flag); + fix[list_initial_integrate_respa[i]]-> + initial_integrate_respa(vflag,ilevel,flag); } /* ---------------------------------------------------------------------- diff --git a/src/modify.h b/src/modify.h index 79f6cbddba..7398a8cef9 100644 --- a/src/modify.h +++ b/src/modify.h @@ -37,8 +37,8 @@ class Modify : protected Pointers { Modify(class LAMMPS *); ~Modify(); void init(); - void setup(); - void initial_integrate(); + void setup(int); + void initial_integrate(int); void pre_decide(); void pre_exchange(); void pre_neighbor(); @@ -47,7 +47,7 @@ class Modify : protected Pointers { void end_of_step(); double thermo_energy(); - void initial_integrate_respa(int,int); + void initial_integrate_respa(int,int,int); void post_force_respa(int,int,int); void final_integrate_respa(int); diff --git a/src/respa.cpp b/src/respa.cpp index e201d22a23..d8df6e9243 100644 --- a/src/respa.cpp +++ b/src/respa.cpp @@ -357,7 +357,7 @@ void Respa::setup() copy_f_flevel(ilevel); } - modify->setup(); + modify->setup(vflag); sum_flevel_f(); output->setup(1); } @@ -373,7 +373,6 @@ void Respa::iterate(int n) for (int i = 0; i < n; i++) { ntimestep = ++update->ntimestep; - ev_set(ntimestep); recurse(nlevels-1); @@ -415,7 +414,7 @@ void Respa::recurse(int ilevel) for (int iloop = 0; iloop < loop[ilevel]; iloop++) { - modify->initial_integrate_respa(ilevel,0); + modify->initial_integrate_respa(vflag,ilevel,0); if (ilevel) recurse(ilevel-1); @@ -427,7 +426,7 @@ void Respa::recurse(int ilevel) // at middle levels, do nothing if (ilevel == nlevels-1) { - modify->initial_integrate_respa(ilevel,1); + modify->initial_integrate_respa(vflag,ilevel,1); int nflag = neighbor->decide(); if (nflag) { if (modify->n_pre_exchange) modify->pre_exchange(); diff --git a/src/verlet.cpp b/src/verlet.cpp index 4a3872fc0e..6e6c66d6f2 100644 --- a/src/verlet.cpp +++ b/src/verlet.cpp @@ -117,7 +117,7 @@ void Verlet::setup() if (force->newton) comm->reverse_communicate(); - modify->setup(); + modify->setup(vflag); output->setup(1); } @@ -132,10 +132,11 @@ void Verlet::iterate(int n) for (int i = 0; i < n; i++) { ntimestep = ++update->ntimestep; + ev_set(ntimestep); // initial time integration - modify->initial_integrate(); + modify->initial_integrate(vflag); // regular communication vs neighbor list rebuild @@ -166,7 +167,6 @@ void Verlet::iterate(int n) // force computations - ev_set(ntimestep); force_clear(); timer->stamp();