From 30f801f5cbc15f9837cd7f20c49e99e0b26418ed Mon Sep 17 00:00:00 2001 From: sjplimp Date: Tue, 8 Apr 2008 22:55:40 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@1725 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/ASPHERE/fix_npt_asphere.cpp | 18 ++-------- src/ASPHERE/fix_npt_asphere.h | 2 -- src/ASPHERE/fix_nve_asphere.cpp | 2 ++ src/ASPHERE/fix_nve_asphere.h | 1 + src/ASPHERE/fix_nvt_asphere.cpp | 18 ++-------- src/ASPHERE/fix_nvt_asphere.h | 2 -- src/atom.cpp | 29 +++++++++------- src/fix_npt.h | 4 +-- src/fix_npt_sphere.cpp | 30 +++++++++------- src/fix_npt_sphere.h | 2 -- src/fix_nve.cpp | 2 -- src/fix_nve.h | 2 +- src/fix_nve_sphere.cpp | 61 +++++++++++++-------------------- src/fix_nve_sphere.h | 10 ++---- src/fix_nvt.h | 2 +- src/fix_nvt_sphere.cpp | 30 +++++++++------- src/fix_nvt_sphere.h | 2 -- src/velocity.cpp | 27 ++++----------- 18 files changed, 93 insertions(+), 151 deletions(-) diff --git a/src/ASPHERE/fix_npt_asphere.cpp b/src/ASPHERE/fix_npt_asphere.cpp index 621c262911..e43de5d95c 100755 --- a/src/ASPHERE/fix_npt_asphere.cpp +++ b/src/ASPHERE/fix_npt_asphere.cpp @@ -44,14 +44,6 @@ FixNPTAsphere::FixNPTAsphere(LAMMPS *lmp, int narg, char **arg) : "quat, angmom, torque, shape"); } -/* ---------------------------------------------------------------------- */ - -void FixNPTAsphere::init() -{ - FixNPT::init(); - dtq = 0.5 * update->dt; -} - /* ---------------------------------------------------------------------- 1st half of Verlet update ------------------------------------------------------------------------- */ @@ -61,6 +53,8 @@ void FixNPTAsphere::initial_integrate(int vflag) int i; double dtfm; + dtq = 0.5 * dtv; + double delta = update->ntimestep - update->beginstep; delta /= update->endstep - update->beginstep; @@ -249,14 +243,6 @@ void FixNPTAsphere::final_integrate() } } -/* ---------------------------------------------------------------------- */ - -void FixNPTAsphere::reset_dt() -{ - FixNPT::reset_dt(); - dtq = 0.5 * update->dt; -} - /* ---------------------------------------------------------------------- Richardson iteration to update quaternion accurately ------------------------------------------------------------------------- */ diff --git a/src/ASPHERE/fix_npt_asphere.h b/src/ASPHERE/fix_npt_asphere.h index 1f6b4420ef..c56e37b7fb 100755 --- a/src/ASPHERE/fix_npt_asphere.h +++ b/src/ASPHERE/fix_npt_asphere.h @@ -21,10 +21,8 @@ namespace LAMMPS_NS { class FixNPTAsphere : public FixNPT { public: FixNPTAsphere(class LAMMPS *, int, char **); - void init(); void initial_integrate(int); void final_integrate(); - void reset_dt(); private: double dtq; diff --git a/src/ASPHERE/fix_nve_asphere.cpp b/src/ASPHERE/fix_nve_asphere.cpp index d3f71e0556..6d13adf675 100755 --- a/src/ASPHERE/fix_nve_asphere.cpp +++ b/src/ASPHERE/fix_nve_asphere.cpp @@ -66,6 +66,8 @@ void FixNVEAsphere::initial_integrate(int vflag) { double dtfm; + dtq = 0.5 * dtv; + double **x = atom->x; double **v = atom->v; double **f = atom->f; diff --git a/src/ASPHERE/fix_nve_asphere.h b/src/ASPHERE/fix_nve_asphere.h index 381084e772..995a29cd14 100755 --- a/src/ASPHERE/fix_nve_asphere.h +++ b/src/ASPHERE/fix_nve_asphere.h @@ -27,6 +27,7 @@ class FixNVEAsphere : public FixNVE { void final_integrate(); private: + double dtq; double **inertia; void richardson(double *, double *, double *); diff --git a/src/ASPHERE/fix_nvt_asphere.cpp b/src/ASPHERE/fix_nvt_asphere.cpp index 0bc410c5aa..eaad9c7544 100755 --- a/src/ASPHERE/fix_nvt_asphere.cpp +++ b/src/ASPHERE/fix_nvt_asphere.cpp @@ -47,18 +47,12 @@ FixNVTAsphere::FixNVTAsphere(LAMMPS *lmp, int narg, char **arg) : /* ---------------------------------------------------------------------- */ -void FixNVTAsphere::init() -{ - FixNVT::init(); - dtq = 0.5 * update->dt; -} - -/* ---------------------------------------------------------------------- */ - void FixNVTAsphere::initial_integrate(int vflag) { double dtfm; + dtq = 0.5 * dtv; + double delta = update->ntimestep - update->beginstep; delta /= update->endstep - update->beginstep; t_target = t_start + delta * (t_stop-t_start); @@ -197,14 +191,6 @@ void FixNVTAsphere::final_integrate() eta_dot *= drag_factor; } -/* ---------------------------------------------------------------------- */ - -void FixNVTAsphere::reset_dt() -{ - FixNVT::reset_dt(); - dtq = 0.5 * update->dt; -} - /* ---------------------------------------------------------------------- Richardson iteration to update quaternion accurately ------------------------------------------------------------------------- */ diff --git a/src/ASPHERE/fix_nvt_asphere.h b/src/ASPHERE/fix_nvt_asphere.h index a53f054a3c..9a108900f3 100755 --- a/src/ASPHERE/fix_nvt_asphere.h +++ b/src/ASPHERE/fix_nvt_asphere.h @@ -21,10 +21,8 @@ namespace LAMMPS_NS { class FixNVTAsphere : public FixNVT { public: FixNVTAsphere(class LAMMPS *, int, char **); - void init(); void initial_integrate(int); void final_integrate(); - void reset_dt(); private: double dtq; diff --git a/src/atom.cpp b/src/atom.cpp index 7f8b7b274d..afd30b6a98 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -601,26 +601,29 @@ void Atom::tag_extend() } /* ---------------------------------------------------------------------- - check if atom tags are consecutive from 1 to Natoms - return 0 if any tag <= 0 or maxtag != Natoms - return 1 if OK (doesn't actually check if all tag values are used) + check that atom IDs span range from 1 to Natoms + return 0 if mintag != 1 or maxtag != Natoms + return 1 if OK + doesn't actually check if all tag values are used ------------------------------------------------------------------------- */ int Atom::tag_consecutive() { - // check[0] = flagged if any tag <= 0 - // check[1] = max tag of any atom - int check[2],check_all[2]; - check[0] = check[1] = 0; - for (int i = 0; i < nlocal; i++) { - if (tag[i] <= 0) check[0] = 1; - if (tag[i] > check[1]) check[1] = tag[i]; - } - MPI_Allreduce(check,check_all,2,MPI_INT,MPI_MAX,world); - if (check_all[0] || check_all[1] != static_cast (natoms)) return 0; + int idmin = static_cast (natoms); + int idmax = 0; + + for (int i = 0; i < nlocal; i++) { + idmin = MIN(idmin,tag[i]); + idmax = MAX(idmax,tag[i]); + } + int idminall,idmaxall; + MPI_Allreduce(&idmin,&idminall,1,MPI_INT,MPI_MIN,world); + MPI_Allreduce(&idmax,&idmaxall,1,MPI_INT,MPI_MAX,world); + + if (idminall != 1 || idmaxall != static_cast (natoms)) return 0; return 1; } diff --git a/src/fix_npt.h b/src/fix_npt.h index c770e9fd16..181c3ed3f9 100644 --- a/src/fix_npt.h +++ b/src/fix_npt.h @@ -25,7 +25,7 @@ class FixNPT : public Fix { int setmask(); virtual void init(); void setup(int); - void initial_integrate(int); + virtual void initial_integrate(int); virtual void final_integrate(); void initial_integrate_respa(int, int, int); void final_integrate_respa(int); @@ -33,7 +33,7 @@ class FixNPT : public Fix { void write_restart(FILE *); void restart(char *); int modify_param(int, char **); - virtual void reset_dt(); + void reset_dt(); protected: int dimension,which; diff --git a/src/fix_npt_sphere.cpp b/src/fix_npt_sphere.cpp index 23894b13d1..bdeb0a5ee7 100644 --- a/src/fix_npt_sphere.cpp +++ b/src/fix_npt_sphere.cpp @@ -53,18 +53,14 @@ FixNPTSphere::~FixNPTSphere() void FixNPTSphere::init() { FixNPT::init(); - dtfrotate = dtf / INERTIA; if (!atom->shape) error->all("Fix npt/sphere requires atom attribute shape"); - double *mass = atom->mass; double **shape = atom->shape; - for (int i = 1; i <= atom->ntypes; i++) { + for (int i = 1; i <= atom->ntypes; i++) if (shape[i][0] != shape[i][1] || shape[i][0] != shape[i][2]) error->all("Fix npt/sphere requires spherical particle shapes"); - dttype[i] = dtfrotate / (0.25*shape[i][0]*shape[i][0]*mass[i]); - } } /* ---------------------------------------------------------------------- @@ -155,6 +151,14 @@ void FixNPTSphere::initial_integrate(int vflag) } } + // recompute timesteps since dt may have changed or come via rRESPA + + double dtfrotate = dtf / INERTIA; + int ntypes = atom->ntypes; + double **shape = atom->shape; + for (int i = 1; i <= ntypes; i++) + dttype[i] = dtfrotate / (0.25*shape[i][0]*shape[i][0]*mass[i]); + // update angular momentum by 1/2 step // update quaternion a full step via Richardson iteration // returns new normalized quaternion @@ -196,6 +200,14 @@ void FixNPTSphere::final_integrate() int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; + // recompute timesteps since dt may have changed or come via rRESPA + + double dtfrotate = dtf / INERTIA; + int ntypes = atom->ntypes; + double **shape = atom->shape; + for (int i = 1; i <= ntypes; i++) + dttype[i] = dtfrotate / (0.25*shape[i][0]*shape[i][0]*mass[i]); + if (which == NOBIAS) { for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { @@ -268,11 +280,3 @@ void FixNPTSphere::final_integrate() omega_dot[i] *= drag_factor; } } - -/* ---------------------------------------------------------------------- */ - -void FixNPTSphere::reset_dt() -{ - FixNPT::reset_dt(); - dtfrotate = dtf / INERTIA; -} diff --git a/src/fix_npt_sphere.h b/src/fix_npt_sphere.h index 4424e14a0f..a31979da8f 100644 --- a/src/fix_npt_sphere.h +++ b/src/fix_npt_sphere.h @@ -25,10 +25,8 @@ class FixNPTSphere : public FixNPT { void init(); void initial_integrate(int); void final_integrate(); - void reset_dt(); private: - double dtfrotate; double *dttype; double factor_rotate; }; diff --git a/src/fix_nve.cpp b/src/fix_nve.cpp index 52030ac7d3..bfe5e74623 100644 --- a/src/fix_nve.cpp +++ b/src/fix_nve.cpp @@ -50,7 +50,6 @@ void FixNVE::init() { dtv = update->dt; dtf = 0.5 * update->dt * force->ftm2v; - dtq = 0.5 * update->dt; if (strcmp(update->integrate_style,"respa") == 0) step_respa = ((Respa *) update->integrate)->step; @@ -169,5 +168,4 @@ void FixNVE::reset_dt() { dtv = update->dt; dtf = 0.5 * update->dt * force->ftm2v; - dtq = 0.5 * update->dt; } diff --git a/src/fix_nve.h b/src/fix_nve.h index 88acbac073..35e7b02f69 100644 --- a/src/fix_nve.h +++ b/src/fix_nve.h @@ -30,7 +30,7 @@ class FixNVE : public Fix { void reset_dt(); protected: - double dtv,dtf,dtq; + double dtv,dtf; double *step_respa; int mass_require; }; diff --git a/src/fix_nve_sphere.cpp b/src/fix_nve_sphere.cpp index f960aef846..a049332e75 100644 --- a/src/fix_nve_sphere.cpp +++ b/src/fix_nve_sphere.cpp @@ -31,7 +31,7 @@ enum{NONE,DIPOLE}; /* ---------------------------------------------------------------------- */ FixNVESphere::FixNVESphere(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg) + FixNVE(lmp, narg, arg) { if (narg < 3) error->all("Illegal fix nve/sphere command"); @@ -86,7 +86,6 @@ void FixNVESphere::init() { dtv = update->dt; dtf = 0.5 * update->dt * force->ftm2v; - dtfrotate = dtf / INERTIA; if (strcmp(update->integrate_style,"respa") == 0) step_respa = ((Respa *) update->integrate)->step; @@ -97,13 +96,10 @@ void FixNVESphere::init() error->all("Fix nve/sphere requires atom attributes radius, rmass"); if (atom->mass) { - double *mass = atom->mass; double **shape = atom->shape; - for (int i = 1; i <= atom->ntypes; i++) { + for (int i = 1; i <= atom->ntypes; i++) if (shape[i][0] != shape[i][1] || shape[i][0] != shape[i][2]) error->all("Fix nve/sphere requires spherical particle shapes"); - dttype[i] = dtfrotate / (0.25*shape[i][0]*shape[i][0]*mass[i]); - } } } @@ -128,6 +124,16 @@ void FixNVESphere::initial_integrate(int vflag) int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; + // recompute timesteps since dt may have changed or come via rRESPA + + dtfrotate = dtf / INERTIA; + if (mass) { + double **shape = atom->shape; + int ntypes = atom->ntypes; + for (int i = 1; i <= ntypes; i++) + dttype[i] = dtfrotate / (0.25*shape[i][0]*shape[i][0]*mass[i]); + } + // update v,x,omega for all particles // d_omega/dt = torque / inertia @@ -212,6 +218,16 @@ void FixNVESphere::final_integrate() int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; + // recompute timesteps since dt may have changed or come via rRESPA + + dtfrotate = dtf / INERTIA; + if (mass) { + double **shape = atom->shape; + int ntypes = atom->ntypes; + for (int i = 1; i <= ntypes; i++) + dttype[i] = dtfrotate / (0.25*shape[i][0]*shape[i][0]*mass[i]); + } + if (mass) { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { @@ -229,6 +245,7 @@ void FixNVESphere::final_integrate() } } else { + dtfrotate = dtf / INERTIA; for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { dtfm = dtf / rmass[i]; @@ -244,35 +261,3 @@ void FixNVESphere::final_integrate() } } } - -/* ---------------------------------------------------------------------- */ - -void FixNVESphere::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; - dtfrotate = dtf / INERTIA; - - if (ilevel == 0) initial_integrate(vflag); - else final_integrate(); -} - -/* ---------------------------------------------------------------------- */ - -void FixNVESphere::final_integrate_respa(int ilevel) -{ - dtf = 0.5 * step_respa[ilevel] * force->ftm2v; - dtfrotate = dtf / INERTIA; - final_integrate(); -} - -/* ---------------------------------------------------------------------- */ - -void FixNVESphere::reset_dt() -{ - dtv = update->dt; - dtf = 0.5 * update->dt * force->ftm2v; - dtfrotate = dtf / INERTIA; -} diff --git a/src/fix_nve_sphere.h b/src/fix_nve_sphere.h index a8b1e8a94c..72b2ddd1e4 100644 --- a/src/fix_nve_sphere.h +++ b/src/fix_nve_sphere.h @@ -14,11 +14,11 @@ #ifndef FIX_NVE_SPHERE_H #define FIX_NVE_SPHERE_H -#include "fix.h" +#include "fix_nve.h" namespace LAMMPS_NS { -class FixNVESphere : public Fix { +class FixNVESphere : public FixNVE { public: FixNVESphere(class LAMMPS *, int, char **); ~FixNVESphere(); @@ -26,14 +26,10 @@ class FixNVESphere : public Fix { void init(); void initial_integrate(int); void final_integrate(); - void initial_integrate_respa(int, int, int); - void final_integrate_respa(int); - void reset_dt(); private: int extra; - double dtv,dtf,dtfrotate; - double *step_respa; + double dtfrotate; double *dttype; }; diff --git a/src/fix_nvt.h b/src/fix_nvt.h index b172813a35..2e0d634966 100644 --- a/src/fix_nvt.h +++ b/src/fix_nvt.h @@ -34,7 +34,7 @@ class FixNVT : public Fix { void restart(char *); int modify_param(int, char **); void reset_target(double); - virtual void reset_dt(); + void reset_dt(); protected: int which; diff --git a/src/fix_nvt_sphere.cpp b/src/fix_nvt_sphere.cpp index 61cadb8157..67faeafb3d 100644 --- a/src/fix_nvt_sphere.cpp +++ b/src/fix_nvt_sphere.cpp @@ -54,18 +54,14 @@ FixNVTSphere::~FixNVTSphere() void FixNVTSphere::init() { FixNVT::init(); - dtfrotate = dtf / INERTIA; if (!atom->shape) error->all("Fix nvt/sphere requires atom attribute shape"); - double *mass = atom->mass; double **shape = atom->shape; - for (int i = 1; i <= atom->ntypes; i++) { + for (int i = 1; i <= atom->ntypes; i++) if (shape[i][0] != shape[i][1] || shape[i][0] != shape[i][2]) error->all("Fix nvt/sphere requires spherical particle shapes"); - dttype[i] = dtfrotate / (0.25*shape[i][0]*shape[i][0]*mass[i]); - } } /* ---------------------------------------------------------------------- */ @@ -101,6 +97,14 @@ void FixNVTSphere::initial_integrate(int vflag) int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; + // recompute timesteps since dt may have changed or come via rRESPA + + double dtfrotate = dtf / INERTIA; + int ntypes = atom->ntypes; + double **shape = atom->shape; + for (int i = 1; i <= ntypes; i++) + dttype[i] = dtfrotate / (0.25*shape[i][0]*shape[i][0]*mass[i]); + if (which == NOBIAS) { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { @@ -164,6 +168,14 @@ void FixNVTSphere::final_integrate() int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; + // recompute timesteps since dt may have changed or come via rRESPA + + double dtfrotate = dtf / INERTIA; + int ntypes = atom->ntypes; + double **shape = atom->shape; + for (int i = 1; i <= ntypes; i++) + dttype[i] = dtfrotate / (0.25*shape[i][0]*shape[i][0]*mass[i]); + if (which == NOBIAS) { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { @@ -211,11 +223,3 @@ void FixNVTSphere::final_integrate() eta_dot += f_eta*dthalf; eta_dot *= drag_factor; } - -/* ---------------------------------------------------------------------- */ - -void FixNVTSphere::reset_dt() -{ - FixNVT::reset_dt(); - dtfrotate = dtf / INERTIA; -} diff --git a/src/fix_nvt_sphere.h b/src/fix_nvt_sphere.h index 23e0e9d91f..6a05f8248b 100644 --- a/src/fix_nvt_sphere.h +++ b/src/fix_nvt_sphere.h @@ -25,10 +25,8 @@ class FixNVTSphere : public FixNVT { void init(); void initial_integrate(int); void final_integrate(); - void reset_dt(); private: - double dtfrotate; double *dttype; }; diff --git a/src/velocity.cpp b/src/velocity.cpp index 7255a7076d..1b10d66fa6 100644 --- a/src/velocity.cpp +++ b/src/velocity.cpp @@ -197,35 +197,20 @@ void Velocity::create(int narg, char **arg) atom->map_set(); } - random = new RanPark(lmp,seed); + // error check if (atom->tag_enable == 0) error->all("Cannot use velocity create loop all unless atoms have IDs"); - int natoms = static_cast (atom->natoms); + if (atom->tag_consecutive() == 0) + error->all("Atom IDs must be consecutive for velocity create loop all"); - // check that atom IDs span range from 1 to natoms - - int *tag = atom->tag; - int idmin = natoms; - int idmax = 0; - - for (i = 0; i < nlocal; i++) { - idmin = MIN(idmin,tag[i]); - idmax = MAX(idmax,tag[i]); - } - int idminall,idmaxall; - MPI_Allreduce(&idmin,&idminall,1,MPI_INT,MPI_MIN,world); - MPI_Allreduce(&idmax,&idmaxall,1,MPI_INT,MPI_MAX,world); - - if (idminall != 1 || idmaxall != natoms) { - char *str = (char *) "Cannot use velocity create loop all with non-contiguous atom IDs"; - error->all(str); - } - // loop over all atoms in system // generate RNGs for all atoms, only assign to ones I own // use either per-type mass or per-atom rmass + random = new RanPark(lmp,seed); + int natoms = static_cast (atom->natoms); + for (i = 1; i <= natoms; i++) { if (dist_flag == 0) { vx = random->uniform();