diff --git a/src/ASPHERE/compute_erotate_asphere.cpp b/src/ASPHERE/compute_erotate_asphere.cpp index 1a6b705856..612619f810 100644 --- a/src/ASPHERE/compute_erotate_asphere.cpp +++ b/src/ASPHERE/compute_erotate_asphere.cpp @@ -25,7 +25,7 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -ComputeERotateASphere::ComputeERotateASphere(LAMMPS *lmp, int narg, char **arg) : +ComputeERotateAsphere::ComputeERotateAsphere(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { if (narg != 3) error->all("Illegal compute erotate/asphere command"); @@ -42,14 +42,14 @@ ComputeERotateASphere::ComputeERotateASphere(LAMMPS *lmp, int narg, char **arg) /* ---------------------------------------------------------------------- */ -ComputeERotateASphere::~ComputeERotateASphere() +ComputeERotateAsphere::~ComputeERotateAsphere() { memory->destroy_2d_double_array(inertia); } /* ---------------------------------------------------------------------- */ -void ComputeERotateASphere::init() +void ComputeERotateAsphere::init() { pfactor = 0.5 * force->mvv2e; @@ -61,7 +61,7 @@ void ComputeERotateASphere::init() /* ---------------------------------------------------------------------- */ -double ComputeERotateASphere::compute_scalar() +double ComputeERotateAsphere::compute_scalar() { invoked |= INVOKED_SCALAR; @@ -102,7 +102,7 @@ double ComputeERotateASphere::compute_scalar() principal moments of inertia for ellipsoids ------------------------------------------------------------------------- */ -void ComputeERotateASphere::calculate_inertia() +void ComputeERotateAsphere::calculate_inertia() { double *mass = atom->mass; double **shape = atom->shape; diff --git a/src/ASPHERE/compute_erotate_asphere.h b/src/ASPHERE/compute_erotate_asphere.h index b098af8519..f058822fb8 100644 --- a/src/ASPHERE/compute_erotate_asphere.h +++ b/src/ASPHERE/compute_erotate_asphere.h @@ -18,10 +18,10 @@ namespace LAMMPS_NS { -class ComputeERotateASphere : public Compute { +class ComputeERotateAsphere : public Compute { public: - ComputeERotateASphere(class LAMMPS *, int, char **); - ~ComputeERotateASphere(); + ComputeERotateAsphere(class LAMMPS *, int, char **); + ~ComputeERotateAsphere(); void init(); double compute_scalar(); diff --git a/src/ASPHERE/compute_temp_asphere.cpp b/src/ASPHERE/compute_temp_asphere.cpp index a8409129d2..e39b82652d 100755 --- a/src/ASPHERE/compute_temp_asphere.cpp +++ b/src/ASPHERE/compute_temp_asphere.cpp @@ -70,6 +70,13 @@ void ComputeTempAsphere::init() fix_dof += modify->fix[i]->dof(igroup); recount(); + if (id_bias) { + tempbias = 1; + int i = modify->find_compute(id_bias); + if (i < 0) error->all("Could not find compute ID for temperature bias"); + tbias = modify->compute[i]; + } else tempbias = 0; + calculate_inertia(); } @@ -122,6 +129,12 @@ double ComputeTempAsphere::compute_scalar() { invoked |= INVOKED_SCALAR; + if (tbias) { + if (!(tbias->invoked & INVOKED_SCALAR)) + double tmp = tbias->compute_scalar(); + tbias->remove_bias_all(); + } + double **v = atom->v; double **quat = atom->quat; double **angmom = atom->angmom; @@ -163,6 +176,8 @@ double ComputeTempAsphere::compute_scalar() } } + if (tbias) tbias->restore_bias_all(); + MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic) recount(); scalar *= tfactor; @@ -177,6 +192,11 @@ void ComputeTempAsphere::compute_vector() invoked |= INVOKED_VECTOR; + if (tbias) { + if (!(tbias->invoked & INVOKED_VECTOR)) tbias->compute_vector(); + tbias->remove_bias_all(); + } + double **v = atom->v; double **quat = atom->quat; double **angmom = atom->angmom; @@ -223,6 +243,8 @@ void ComputeTempAsphere::compute_vector() t[5] += inertia[itype][2]*wbody[1]*wbody[2]; } + if (tbias) tbias->restore_bias_all(); + MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); for (i = 0; i < 6; i++) vector[i] *= force->mvv2e; } @@ -245,3 +267,42 @@ void ComputeTempAsphere::calculate_inertia() (shape[i][0]*shape[i][0]+shape[i][1]*shape[i][1]) / 5.0; } } + +/* ---------------------------------------------------------------------- + remove velocity bias from atom I to leave thermal velocity +------------------------------------------------------------------------- */ + +void ComputeTempAsphere::remove_bias(int i, double *v) +{ + if (tbias) tbias->remove_bias(i,v); +} + +/* ---------------------------------------------------------------------- + remove velocity bias from all atoms to leave thermal velocity +------------------------------------------------------------------------- */ + +void ComputeTempAsphere::remove_bias_all() +{ + if (tbias) tbias->remove_bias_all(); +} + +/* ---------------------------------------------------------------------- + add back in velocity bias to atom I removed by remove_bias() + assume remove_bias() was previously called +------------------------------------------------------------------------- */ + +void ComputeTempAsphere::restore_bias(double *v) +{ + if (tbias) tbias->restore_bias(v); +} + +/* ---------------------------------------------------------------------- + add back in velocity bias to all atoms removed by remove_bias_all() + assume remove_bias_all() was previously called +------------------------------------------------------------------------- */ + +void ComputeTempAsphere::restore_bias_all() +{ + if (tbias) tbias->restore_bias_all(); +} + diff --git a/src/ASPHERE/compute_temp_asphere.h b/src/ASPHERE/compute_temp_asphere.h index 651f409853..42441face1 100755 --- a/src/ASPHERE/compute_temp_asphere.h +++ b/src/ASPHERE/compute_temp_asphere.h @@ -26,11 +26,18 @@ class ComputeTempAsphere : public Compute { double compute_scalar(); void compute_vector(); + void remove_bias(int, double *); + void remove_bias_all(); + void restore_bias(double *); + void restore_bias_all(); + private: int fix_dof; double tfactor; double **inertia; + Compute *tbias; // ptr to additional bias compute + void recount(); void calculate_inertia(); }; diff --git a/src/ASPHERE/fix_npt_asphere.cpp b/src/ASPHERE/fix_npt_asphere.cpp index a5f2752bcc..edb2897378 100755 --- a/src/ASPHERE/fix_npt_asphere.cpp +++ b/src/ASPHERE/fix_npt_asphere.cpp @@ -35,7 +35,7 @@ enum{NOBIAS,BIAS}; /* ---------------------------------------------------------------------- */ -FixNPTASphere::FixNPTASphere(LAMMPS *lmp, int narg, char **arg) : +FixNPTAsphere::FixNPTAsphere(LAMMPS *lmp, int narg, char **arg) : FixNPT(lmp, narg, arg) { if (!atom->quat_flag || !atom->angmom_flag || !atom->torque_flag || @@ -46,7 +46,7 @@ FixNPTASphere::FixNPTASphere(LAMMPS *lmp, int narg, char **arg) : /* ---------------------------------------------------------------------- */ -void FixNPTASphere::init() +void FixNPTAsphere::init() { FixNPT::init(); dtq = 0.5 * update->dt; @@ -56,7 +56,7 @@ void FixNPTASphere::init() 1st half of Verlet update ------------------------------------------------------------------------- */ -void FixNPTASphere::initial_integrate(int vflag) +void FixNPTAsphere::initial_integrate(int vflag) { int i; double dtfm; @@ -168,7 +168,7 @@ void FixNPTASphere::initial_integrate(int vflag) 2nd half of Verlet update ------------------------------------------------------------------------- */ -void FixNPTASphere::final_integrate() +void FixNPTAsphere::final_integrate() { int i; double dtfm; @@ -251,7 +251,7 @@ void FixNPTASphere::final_integrate() /* ---------------------------------------------------------------------- */ -void FixNPTASphere::reset_dt() +void FixNPTAsphere::reset_dt() { FixNPT::reset_dt(); dtq = 0.5 * update->dt; @@ -261,7 +261,7 @@ void FixNPTASphere::reset_dt() Richardson iteration to update quaternion accurately ------------------------------------------------------------------------- */ -void FixNPTASphere::richardson(double *q, double *m, double *moments) +void FixNPTAsphere::richardson(double *q, double *m, double *moments) { // compute omega at 1/2 step from m at 1/2 step and q at 0 @@ -320,7 +320,7 @@ void FixNPTASphere::richardson(double *q, double *m, double *moments) and divide by principal moments ------------------------------------------------------------------------- */ -void FixNPTASphere::omega_from_mq(double *q, double *m, double *inertia, +void FixNPTAsphere::omega_from_mq(double *q, double *m, double *inertia, double *w) { double rot[3][3]; @@ -339,7 +339,7 @@ void FixNPTASphere::omega_from_mq(double *q, double *m, double *inertia, shape = x,y,z radii in body frame ------------------------------------------------------------------------- */ -void FixNPTASphere::calculate_inertia(double mass, double *shape, +void FixNPTAsphere::calculate_inertia(double mass, double *shape, double *inertia) { inertia[0] = mass*(shape[1]*shape[1]+shape[2]*shape[2])/5.0; diff --git a/src/ASPHERE/fix_npt_asphere.h b/src/ASPHERE/fix_npt_asphere.h index e7e19b07f4..1f6b4420ef 100755 --- a/src/ASPHERE/fix_npt_asphere.h +++ b/src/ASPHERE/fix_npt_asphere.h @@ -18,9 +18,9 @@ namespace LAMMPS_NS { -class FixNPTASphere : public FixNPT { +class FixNPTAsphere : public FixNPT { public: - FixNPTASphere(class LAMMPS *, int, char **); + FixNPTAsphere(class LAMMPS *, int, char **); void init(); void initial_integrate(int); void final_integrate(); diff --git a/src/ASPHERE/fix_nve_asphere.cpp b/src/ASPHERE/fix_nve_asphere.cpp index 5286913229..d3f71e0556 100755 --- a/src/ASPHERE/fix_nve_asphere.cpp +++ b/src/ASPHERE/fix_nve_asphere.cpp @@ -34,7 +34,7 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -FixNVEASphere::FixNVEASphere(LAMMPS *lmp, int narg, char **arg) : +FixNVEAsphere::FixNVEAsphere(LAMMPS *lmp, int narg, char **arg) : FixNVE(lmp, narg, arg) { if (!atom->quat_flag || !atom->angmom_flag || !atom->torque_flag || @@ -47,14 +47,14 @@ FixNVEASphere::FixNVEASphere(LAMMPS *lmp, int narg, char **arg) : /* ---------------------------------------------------------------------- */ -FixNVEASphere::~FixNVEASphere() +FixNVEAsphere::~FixNVEAsphere() { memory->destroy_2d_double_array(inertia); } /* ---------------------------------------------------------------------- */ -void FixNVEASphere::init() +void FixNVEAsphere::init() { FixNVE::init(); calculate_inertia(); @@ -62,7 +62,7 @@ void FixNVEASphere::init() /* ---------------------------------------------------------------------- */ -void FixNVEASphere::initial_integrate(int vflag) +void FixNVEAsphere::initial_integrate(int vflag) { double dtfm; @@ -102,7 +102,7 @@ void FixNVEASphere::initial_integrate(int vflag) /* ---------------------------------------------------------------------- */ -void FixNVEASphere::final_integrate() +void FixNVEAsphere::final_integrate() { double dtfm; @@ -133,7 +133,7 @@ void FixNVEASphere::final_integrate() Richardson iteration to update quaternion accurately ------------------------------------------------------------------------- */ -void FixNVEASphere::richardson(double *q, double *m, double *moments) +void FixNVEAsphere::richardson(double *q, double *m, double *moments) { // compute omega at 1/2 step from m at 1/2 step and q at 0 @@ -192,7 +192,7 @@ void FixNVEASphere::richardson(double *q, double *m, double *moments) and divide by principal moments ------------------------------------------------------------------------- */ -void FixNVEASphere::omega_from_mq(double *q, double *m, double *moments, +void FixNVEAsphere::omega_from_mq(double *q, double *m, double *moments, double *w) { double rot[3][3]; @@ -210,7 +210,7 @@ void FixNVEASphere::omega_from_mq(double *q, double *m, double *moments, principal moments of inertia for ellipsoids ------------------------------------------------------------------------- */ -void FixNVEASphere::calculate_inertia() +void FixNVEAsphere::calculate_inertia() { double *mass = atom->mass; double **shape = atom->shape; diff --git a/src/ASPHERE/fix_nve_asphere.h b/src/ASPHERE/fix_nve_asphere.h index c748b6db5e..381084e772 100755 --- a/src/ASPHERE/fix_nve_asphere.h +++ b/src/ASPHERE/fix_nve_asphere.h @@ -18,10 +18,10 @@ namespace LAMMPS_NS { -class FixNVEASphere : public FixNVE { +class FixNVEAsphere : public FixNVE { public: - FixNVEASphere(class LAMMPS *, int, char **); - ~FixNVEASphere(); + FixNVEAsphere(class LAMMPS *, int, char **); + ~FixNVEAsphere(); void init(); void initial_integrate(int); void final_integrate(); diff --git a/src/ASPHERE/fix_nvt_asphere.cpp b/src/ASPHERE/fix_nvt_asphere.cpp index 2b29a7df0a..1c820c766d 100755 --- a/src/ASPHERE/fix_nvt_asphere.cpp +++ b/src/ASPHERE/fix_nvt_asphere.cpp @@ -36,7 +36,7 @@ enum{NOBIAS,BIAS}; /* ---------------------------------------------------------------------- */ -FixNVTASphere::FixNVTASphere(LAMMPS *lmp, int narg, char **arg) : +FixNVTAsphere::FixNVTAsphere(LAMMPS *lmp, int narg, char **arg) : FixNVT(lmp, narg, arg) { if (!atom->quat_flag || !atom->angmom_flag || !atom->torque_flag || @@ -47,7 +47,7 @@ FixNVTASphere::FixNVTASphere(LAMMPS *lmp, int narg, char **arg) : /* ---------------------------------------------------------------------- */ -void FixNVTASphere::init() +void FixNVTAsphere::init() { FixNVT::init(); dtq = 0.5 * update->dt; @@ -55,7 +55,7 @@ void FixNVTASphere::init() /* ---------------------------------------------------------------------- */ -void FixNVTASphere::initial_integrate(int vflag) +void FixNVTAsphere::initial_integrate(int vflag) { double dtfm; @@ -141,7 +141,7 @@ void FixNVTASphere::initial_integrate(int vflag) /* ---------------------------------------------------------------------- */ -void FixNVTASphere::final_integrate() +void FixNVTAsphere::final_integrate() { double dtfm; @@ -199,7 +199,7 @@ void FixNVTASphere::final_integrate() /* ---------------------------------------------------------------------- */ -void FixNVTASphere::reset_dt() +void FixNVTAsphere::reset_dt() { FixNVT::reset_dt(); dtq = 0.5 * update->dt; @@ -209,7 +209,7 @@ void FixNVTASphere::reset_dt() Richardson iteration to update quaternion accurately ------------------------------------------------------------------------- */ -void FixNVTASphere::richardson(double *q, double *m, double *moments) +void FixNVTAsphere::richardson(double *q, double *m, double *moments) { // compute omega at 1/2 step from m at 1/2 step and q at 0 @@ -268,7 +268,7 @@ void FixNVTASphere::richardson(double *q, double *m, double *moments) and divide by principal moments ------------------------------------------------------------------------- */ -void FixNVTASphere::omega_from_mq(double *q, double *m, double *inertia, +void FixNVTAsphere::omega_from_mq(double *q, double *m, double *inertia, double *w) { double rot[3][3]; @@ -287,7 +287,7 @@ void FixNVTASphere::omega_from_mq(double *q, double *m, double *inertia, shape = x,y,z radii in body frame ------------------------------------------------------------------------- */ -void FixNVTASphere::calculate_inertia(double mass, double *shape, +void FixNVTAsphere::calculate_inertia(double mass, double *shape, double *inertia) { inertia[0] = mass*(shape[1]*shape[1]+shape[2]*shape[2])/5.0; diff --git a/src/ASPHERE/fix_nvt_asphere.h b/src/ASPHERE/fix_nvt_asphere.h index 5dd4de37b5..a53f054a3c 100755 --- a/src/ASPHERE/fix_nvt_asphere.h +++ b/src/ASPHERE/fix_nvt_asphere.h @@ -18,9 +18,9 @@ namespace LAMMPS_NS { -class FixNVTASphere : public FixNVT { +class FixNVTAsphere : public FixNVT { public: - FixNVTASphere(class LAMMPS *, int, char **); + FixNVTAsphere(class LAMMPS *, int, char **); void init(); void initial_integrate(int); void final_integrate(); diff --git a/src/ASPHERE/style_asphere.h b/src/ASPHERE/style_asphere.h index 1a84677e10..70dba98aa4 100644 --- a/src/ASPHERE/style_asphere.h +++ b/src/ASPHERE/style_asphere.h @@ -25,7 +25,7 @@ AtomStyle(ellipsoid,AtomVecEllipsoid) #endif #ifdef ComputeClass -ComputeStyle(erotate/asphere,ComputeERotateASphere) +ComputeStyle(erotate/asphere,ComputeERotateAsphere) ComputeStyle(temp/asphere,ComputeTempAsphere) #endif @@ -36,9 +36,9 @@ ComputeStyle(temp/asphere,ComputeTempAsphere) #endif #ifdef FixClass -FixStyle(nve/asphere,FixNVEASphere) -FixStyle(nvt/asphere,FixNVTASphere) -FixStyle(npt/asphere,FixNPTASphere) +FixStyle(nve/asphere,FixNVEAsphere) +FixStyle(nvt/asphere,FixNVTAsphere) +FixStyle(npt/asphere,FixNPTAsphere) #endif #ifdef PairInclude diff --git a/src/compute.cpp b/src/compute.cpp index 83edb2d933..828e07e7f2 100644 --- a/src/compute.cpp +++ b/src/compute.cpp @@ -60,7 +60,10 @@ Compute::Compute(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) scalar_flag = vector_flag = peratom_flag = 0; tempflag = pressflag = peflag = 0; pressatomflag = peatomflag = 0; + tempbias = 0; + id_bias = NULL; + id_pre = NULL; timeflag = 0; invoked = 0; @@ -83,6 +86,7 @@ Compute::~Compute() { delete [] id; delete [] style; + delete [] id_bias; delete [] id_pre; memory->sfree(tlist); @@ -112,6 +116,16 @@ void Compute::modify_params(int narg, char **arg) else if (strcmp(arg[iarg+1],"yes") == 0) thermoflag = 1; else error->all("Illegal compute_modify command"); iarg += 2; + } else if (strcmp(arg[iarg],"bias") == 0) { + if (iarg+2 > narg) error->all("Illegal compute_modify command"); + delete [] id_bias; + if (strcmp(arg[iarg+1],"NULL") == 0) id_bias = NULL; + else { + int n = strlen(arg[iarg+1]) + 1; + id_bias = new char[n]; + strcpy(id_bias,arg[iarg+1]); + } + iarg += 2; } else error->all("Illegal compute_modify command"); } } @@ -163,15 +177,3 @@ int Compute::matchstep(int ntimestep) } return 0; } - -/* ---------------------------------------------------------------------- - add back in velocity bias removed by remove_bias() to leave thermal temp - assume remove_bias() was previously called for this atom -------------------------------------------------------------------------- */ - -void Compute::restore_bias(double *v) -{ - v[0] += vbias[0]; - v[1] += vbias[1]; - v[2] += vbias[2]; -} diff --git a/src/compute.h b/src/compute.h index cde54e3174..9f887bd100 100644 --- a/src/compute.h +++ b/src/compute.h @@ -45,8 +45,9 @@ class Compute : protected Pointers { int peflag; // 1 if Compute calculates PE (uses Force energies) int peatomflag; // 1 if Compute calculates per-atom PE - int tempbias; // 0/1 if has bias routines for isolating thermal temp - double vbias[3]; // stored velocity bias for one atom + int tempbias; // 0/1 if Compute temp includes + // self or extra bias via compute_modify + char *id_bias; // ID of extra Compute temp that adds bias char *id_pre; // ID of pre-compute the Compute may store @@ -77,7 +78,9 @@ class Compute : protected Pointers { virtual void unpack_reverse_comm(int, int *, double *) {} virtual void remove_bias(int, double *) {} - void restore_bias(double *); + virtual void remove_bias_all() {} + virtual void restore_bias(double *) {} + virtual void restore_bias_all() {} void addstep(int); int matchstep(int); diff --git a/src/compute_temp.cpp b/src/compute_temp.cpp index 7fecab3a4a..8b2a8931c2 100644 --- a/src/compute_temp.cpp +++ b/src/compute_temp.cpp @@ -57,6 +57,13 @@ void ComputeTemp::init() for (int i = 0; i < modify->nfix; i++) fix_dof += modify->fix[i]->dof(igroup); recount(); + + if (id_bias) { + tempbias = 1; + int i = modify->find_compute(id_bias); + if (i < 0) error->all("Could not find compute ID for temperature bias"); + tbias = modify->compute[i]; + } } /* ---------------------------------------------------------------------- */ @@ -76,6 +83,12 @@ double ComputeTemp::compute_scalar() { invoked |= INVOKED_SCALAR; + if (tbias) { + if (!(tbias->invoked & INVOKED_SCALAR)) + double tmp = tbias->compute_scalar(); + tbias->remove_bias_all(); + } else tempbias = 0; + double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; @@ -96,6 +109,8 @@ double ComputeTemp::compute_scalar() t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i]; } + if (tbias) tbias->restore_bias_all(); + MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic) recount(); scalar *= tfactor; @@ -110,6 +125,11 @@ void ComputeTemp::compute_vector() invoked |= INVOKED_VECTOR; + if (tbias) { + if (!(tbias->invoked & INVOKED_VECTOR)) tbias->compute_vector(); + tbias->remove_bias_all(); + } + double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; @@ -132,6 +152,46 @@ void ComputeTemp::compute_vector() t[5] += massone * v[i][1]*v[i][2]; } + if (tbias) tbias->restore_bias_all(); + MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); for (i = 0; i < 6; i++) vector[i] *= force->mvv2e; } + +/* ---------------------------------------------------------------------- + remove velocity bias from atom I to leave thermal velocity +------------------------------------------------------------------------- */ + +void ComputeTemp::remove_bias(int i, double *v) +{ + if (tbias) tbias->remove_bias(i,v); +} + +/* ---------------------------------------------------------------------- + remove velocity bias from all atoms to leave thermal velocity +------------------------------------------------------------------------- */ + +void ComputeTemp::remove_bias_all() +{ + if (tbias) tbias->remove_bias_all(); +} + +/* ---------------------------------------------------------------------- + add back in velocity bias to atom I removed by remove_bias() + assume remove_bias() was previously called +------------------------------------------------------------------------- */ + +void ComputeTemp::restore_bias(double *v) +{ + if (tbias) tbias->restore_bias(v); +} + +/* ---------------------------------------------------------------------- + add back in velocity bias to all atoms removed by remove_bias_all() + assume remove_bias_all() was previously called +------------------------------------------------------------------------- */ + +void ComputeTemp::restore_bias_all() +{ + if (tbias) tbias->restore_bias_all(); +} diff --git a/src/compute_temp.h b/src/compute_temp.h index 89db81c534..2405abbd66 100644 --- a/src/compute_temp.h +++ b/src/compute_temp.h @@ -26,10 +26,17 @@ class ComputeTemp : public Compute { double compute_scalar(); void compute_vector(); + void remove_bias(int, double *); + void remove_bias_all(); + void restore_bias(double *); + void restore_bias_all(); + private: int fix_dof; double tfactor; + Compute *tbias; // ptr to additional bias compute + void recount(); }; diff --git a/src/compute_temp_com.cpp b/src/compute_temp_com.cpp index df80cd3873..3b808c7cf2 100644 --- a/src/compute_temp_com.cpp +++ b/src/compute_temp_com.cpp @@ -45,7 +45,6 @@ ComputeTempCOM::ComputeTempCOM(LAMMPS *lmp, int narg, char **arg) : extvector = 1; tempflag = 1; tempbias = 1; - vbias[0] = vbias[1] = vbias[2] = 0.0; vector = new double[6]; } @@ -66,6 +65,12 @@ void ComputeTempCOM::init() fix_dof += modify->fix[i]->dof(igroup); recount(); masstotal = group->mass(igroup); + + if (id_bias) { + int i = modify->find_compute(id_bias); + if (i < 0) error->all("Could not find compute ID for temperature bias"); + tbias = modify->compute[i]; + } } /* ---------------------------------------------------------------------- */ @@ -87,6 +92,12 @@ double ComputeTempCOM::compute_scalar() invoked |= INVOKED_SCALAR; + if (tbias) { + if (!(tbias->invoked & INVOKED_SCALAR)) + double tmp = tbias->compute_scalar(); + tbias->remove_bias_all(); + } + if (dynamic) masstotal = group->mass(igroup); group->vcm(igroup,masstotal,vbias); @@ -111,6 +122,8 @@ double ComputeTempCOM::compute_scalar() vthermal[2]*vthermal[2]) * rmass[i]; } + if (tbias) tbias->restore_bias_all(); + MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic) recount(); scalar *= tfactor; @@ -126,6 +139,11 @@ void ComputeTempCOM::compute_vector() invoked |= INVOKED_VECTOR; + if (tbias) { + if (!(tbias->invoked & INVOKED_VECTOR)) tbias->compute_vector(); + tbias->remove_bias_all(); + } + if (dynamic) masstotal = group->mass(igroup); group->vcm(igroup,masstotal,vbias); @@ -156,15 +174,74 @@ void ComputeTempCOM::compute_vector() t[5] += massone * vthermal[1]*vthermal[2]; } + if (tbias) tbias->restore_bias_all(); + MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); for (i = 0; i < 6; i++) vector[i] *= force->mvv2e; } -/* ---------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + remove velocity bias from atom I to leave thermal velocity +------------------------------------------------------------------------- */ void ComputeTempCOM::remove_bias(int i, double *v) { - v[0] = v[0] - vbias[0]; - v[1] = v[1] - vbias[1]; - v[2] = v[2] - vbias[2]; + if (tbias) tbias->remove_bias(i,v); + v[0] -= vbias[0]; + v[1] -= vbias[1]; + v[2] -= vbias[2]; +} + +/* ---------------------------------------------------------------------- + remove velocity bias from all atoms to leave thermal velocity +------------------------------------------------------------------------- */ + +void ComputeTempCOM::remove_bias_all() +{ + if (tbias) tbias->remove_bias_all(); + + double **v = atom->v; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + v[i][0] -= vbias[0]; + v[i][1] -= vbias[1]; + v[i][2] -= vbias[2]; + } +} + +/* ---------------------------------------------------------------------- + add back in velocity bias to atom I removed by remove_bias() + assume remove_bias() was previously called +------------------------------------------------------------------------- */ + +void ComputeTempCOM::restore_bias(double *v) +{ + v[0] += vbias[0]; + v[1] += vbias[1]; + v[2] += vbias[2]; + if (tbias) tbias->restore_bias(v); +} + +/* ---------------------------------------------------------------------- + add back in velocity bias to all atoms removed by remove_bias_all() + assume remove_bias_all() was previously called +------------------------------------------------------------------------- */ + +void ComputeTempCOM::restore_bias_all() +{ + double **v = atom->v; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + v[i][0] += vbias[0]; + v[i][1] += vbias[1]; + v[i][2] += vbias[2]; + } + + if (tbias) tbias->restore_bias_all(); } diff --git a/src/compute_temp_com.h b/src/compute_temp_com.h index cdb2b261cf..4944b2986c 100644 --- a/src/compute_temp_com.h +++ b/src/compute_temp_com.h @@ -25,12 +25,19 @@ class ComputeTempCOM : public Compute { void init(); double compute_scalar(); void compute_vector(); + void remove_bias(int, double *); + void remove_bias_all(); + void restore_bias(double *); + void restore_bias_all(); private: int fix_dof; double tfactor,masstotal; + double vbias[3]; // stored velocity bias for one atom + Compute *tbias; // ptr to additional bias compute + void recount(); }; diff --git a/src/compute_temp_deform.cpp b/src/compute_temp_deform.cpp index 0413c1a382..cca5e0e5b5 100644 --- a/src/compute_temp_deform.cpp +++ b/src/compute_temp_deform.cpp @@ -26,6 +26,7 @@ #include "fix_deform.h" #include "group.h" #include "comm.h" +#include "memory.h" #include "error.h" using namespace LAMMPS_NS; @@ -48,6 +49,8 @@ ComputeTempDeform::ComputeTempDeform(LAMMPS *lmp, int narg, char **arg) : tempflag = 1; tempbias = 1; + maxbias = 0; + vbiasall = NULL; vector = new double[6]; } @@ -55,6 +58,7 @@ ComputeTempDeform::ComputeTempDeform(LAMMPS *lmp, int narg, char **arg) : ComputeTempDeform::~ComputeTempDeform() { + memory->destroy_2d_double_array(vbiasall); delete [] vector; } @@ -69,6 +73,12 @@ void ComputeTempDeform::init() fix_dof += modify->fix[i]->dof(igroup); recount(); + if (id_bias) { + int i = modify->find_compute(id_bias); + if (i < 0) error->all("Could not find compute ID for temperature bias"); + tbias = modify->compute[i]; + } + // check fix deform remap settings for (i = 0; i < modify->nfix; i++) @@ -101,6 +111,12 @@ double ComputeTempDeform::compute_scalar() invoked |= INVOKED_SCALAR; + if (tbias) { + if (!(tbias->invoked & INVOKED_SCALAR)) + double tmp = tbias->compute_scalar(); + tbias->remove_bias_all(); + } + double **x = atom->x; double **v = atom->v; double *mass = atom->mass; @@ -118,35 +134,25 @@ double ComputeTempDeform::compute_scalar() double t = 0.0; - if (mass) { - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - domain->x2lamda(x[i],lamda); - vstream[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + - h_rate[4]*lamda[2] + h_ratelo[0]; - vstream[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1]; - vstream[2] = h_rate[2]*lamda[2] + h_ratelo[2]; - vthermal[0] = v[i][0] - vstream[0]; - vthermal[1] = v[i][1] - vstream[1]; - vthermal[2] = v[i][2] - vstream[2]; + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + domain->x2lamda(x[i],lamda); + vstream[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + + h_rate[4]*lamda[2] + h_ratelo[0]; + vstream[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1]; + vstream[2] = h_rate[2]*lamda[2] + h_ratelo[2]; + vthermal[0] = v[i][0] - vstream[0]; + vthermal[1] = v[i][1] - vstream[1]; + vthermal[2] = v[i][2] - vstream[2]; + if (mass) t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + vthermal[2]*vthermal[2]) * mass[type[i]]; - } - } - else - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - domain->x2lamda(x[i],lamda); - vstream[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + - h_rate[4]*lamda[2] + h_ratelo[0]; - vstream[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1]; - vstream[2] = h_rate[2]*lamda[2] + h_ratelo[2]; - vthermal[0] = v[i][0] - vstream[0]; - vthermal[1] = v[i][1] - vstream[1]; - vthermal[2] = v[i][2] - vstream[2]; + else t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + vthermal[2]*vthermal[2]) * rmass[i]; - } + } + + if (tbias) tbias->restore_bias_all(); MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); scalar *= tfactor; @@ -161,6 +167,11 @@ void ComputeTempDeform::compute_vector() invoked |= INVOKED_VECTOR; + if (tbias) { + if (!(tbias->invoked & INVOKED_VECTOR)) tbias->compute_vector(); + tbias->remove_bias_all(); + } + double **x = atom->x; double **v = atom->v; double *mass = atom->mass; @@ -196,14 +207,20 @@ void ComputeTempDeform::compute_vector() t[5] += massone * vthermal[1]*vthermal[2]; } + if (tbias) tbias->restore_bias_all(); + MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); for (int i = 0; i < 6; i++) vector[i] *= force->mvv2e; } -/* ---------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + remove velocity bias from atom I to leave thermal velocity +------------------------------------------------------------------------- */ void ComputeTempDeform::remove_bias(int i, double *v) { + if (tbias) tbias->remove_bias(i,v); + double lamda[3]; double *h_rate = domain->h_rate; double *h_ratelo = domain->h_ratelo; @@ -217,3 +234,81 @@ void ComputeTempDeform::remove_bias(int i, double *v) v[1] -= vbias[1]; v[2] -= vbias[2]; } + +/* ---------------------------------------------------------------------- + remove velocity bias from all atoms to leave thermal velocity +------------------------------------------------------------------------- */ + +void ComputeTempDeform::remove_bias_all() +{ + if (tbias) tbias->remove_bias_all(); + + double **v = atom->v; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + if (nlocal > maxbias) { + memory->destroy_2d_double_array(vbiasall); + maxbias = atom->nmax; + vbiasall = memory->create_2d_double_array(maxbias,3, + "compute/temp:vbiasall"); + } + + double lamda[3]; + double *h_rate = domain->h_rate; + double *h_ratelo = domain->h_ratelo; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + domain->x2lamda(atom->x[i],lamda); + vbiasall[i][0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + + h_rate[4]*lamda[2] + h_ratelo[0]; + vbiasall[i][1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1]; + vbiasall[i][2] = h_rate[2]*lamda[2] + h_ratelo[2]; + v[i][0] -= vbiasall[i][0]; + v[i][1] -= vbiasall[i][1]; + v[i][2] -= vbiasall[i][2]; + } +} + +/* ---------------------------------------------------------------------- + add back in velocity bias to atom I removed by remove_bias() + assume remove_bias() was previously called +------------------------------------------------------------------------- */ + +void ComputeTempDeform::restore_bias(double *v) +{ + v[0] += vbias[0]; + v[1] += vbias[1]; + v[2] += vbias[2]; + if (tbias) tbias->restore_bias(v); +} + +/* ---------------------------------------------------------------------- + add back in velocity bias to all atoms removed by remove_bias_all() + assume remove_bias_all() was previously called +------------------------------------------------------------------------- */ + +void ComputeTempDeform::restore_bias_all() +{ + double **v = atom->v; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + v[i][0] += vbiasall[i][0]; + v[i][1] += vbiasall[i][1]; + v[i][2] += vbiasall[i][2]; + } + + if (tbias) tbias->restore_bias_all(); +} + +/* ---------------------------------------------------------------------- */ + +double ComputeTempDeform::memory_usage() +{ + double bytes = maxbias * sizeof(double); + return bytes; +} diff --git a/src/compute_temp_deform.h b/src/compute_temp_deform.h index a489e82765..f6ab938322 100644 --- a/src/compute_temp_deform.h +++ b/src/compute_temp_deform.h @@ -25,12 +25,22 @@ class ComputeTempDeform : public Compute { void init(); double compute_scalar(); void compute_vector(); + void remove_bias(int, double *); + void remove_bias_all(); + void restore_bias(double *); + void restore_bias_all(); + double memory_usage(); private: int fix_dof; double tfactor; + double vbias[3]; // stored velocity bias for one atom + double **vbiasall; // stored velocity bias for all atoms + int maxbias; // size of vbiasall array + Compute *tbias; // ptr to additional bias compute + void recount(); }; diff --git a/src/compute_temp_partial.cpp b/src/compute_temp_partial.cpp index 29034a1d8f..eaacee47e4 100644 --- a/src/compute_temp_partial.cpp +++ b/src/compute_temp_partial.cpp @@ -19,6 +19,7 @@ #include "modify.h" #include "fix.h" #include "group.h" +#include "memory.h" #include "error.h" using namespace LAMMPS_NS; @@ -43,8 +44,9 @@ ComputeTempPartial::ComputeTempPartial(LAMMPS *lmp, int narg, char **arg) : extvector = 1; tempflag = 1; tempbias = 1; - vbias[0] = vbias[1] = vbias[2] = 0.0; + maxbias = 0; + vbiasall = NULL; vector = new double[6]; } @@ -52,6 +54,7 @@ ComputeTempPartial::ComputeTempPartial(LAMMPS *lmp, int narg, char **arg) : ComputeTempPartial::~ComputeTempPartial() { + memory->destroy_2d_double_array(vbiasall); delete [] vector; } @@ -63,6 +66,12 @@ void ComputeTempPartial::init() for (int i = 0; i < modify->nfix; i++) fix_dof += modify->fix[i]->dof(igroup); recount(); + + if (id_bias) { + int i = modify->find_compute(id_bias); + if (i < 0) error->all("Could not find compute ID for temperature bias"); + tbias = modify->compute[i]; + } } /* ---------------------------------------------------------------------- */ @@ -82,6 +91,12 @@ double ComputeTempPartial::compute_scalar() { invoked |= INVOKED_SCALAR; + if (tbias) { + if (!(tbias->invoked & INVOKED_SCALAR)) + double tmp = tbias->compute_scalar(); + tbias->remove_bias_all(); + } + double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; @@ -102,7 +117,9 @@ double ComputeTempPartial::compute_scalar() t += (xflag*v[i][0]*v[i][0] + yflag*v[i][1]*v[i][1] + zflag*v[i][2]*v[i][2]) * rmass[i]; } - + + if (tbias) tbias->restore_bias_all(); + MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic) recount(); scalar *= tfactor; @@ -117,6 +134,11 @@ void ComputeTempPartial::compute_vector() invoked |= INVOKED_VECTOR; + if (tbias) { + if (!(tbias->invoked & INVOKED_VECTOR)) tbias->compute_vector(); + tbias->remove_bias_all(); + } + double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; @@ -139,14 +161,20 @@ void ComputeTempPartial::compute_vector() t[5] += massone * yflag*zflag*v[i][1]*v[i][2]; } + if (tbias) tbias->restore_bias_all(); + MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); for (i = 0; i < 6; i++) vector[i] *= force->mvv2e; } -/* ---------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + remove velocity bias from atom I to leave thermal velocity +------------------------------------------------------------------------- */ void ComputeTempPartial::remove_bias(int i, double *v) { + if (tbias) tbias->remove_bias(i,v); + if (!xflag) { vbias[0] = v[0]; v[0] = 0.0; @@ -160,3 +188,81 @@ void ComputeTempPartial::remove_bias(int i, double *v) v[2] = 0.0; } } + +/* ---------------------------------------------------------------------- + remove velocity bias from all atoms to leave thermal velocity +------------------------------------------------------------------------- */ + +void ComputeTempPartial::remove_bias_all() +{ + if (tbias) tbias->remove_bias_all(); + + double **v = atom->v; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + if (nlocal > maxbias) { + memory->destroy_2d_double_array(vbiasall); + maxbias = atom->nmax; + vbiasall = memory->create_2d_double_array(maxbias,3, + "compute/temp:vbiasall"); + } + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + if (!xflag) { + vbiasall[i][0] = v[i][0]; + v[i][0] = 0.0; + } + if (!yflag) { + vbiasall[i][1] = v[i][1]; + v[i][1] = 0.0; + } + if (!zflag) { + vbiasall[i][2] = v[i][2]; + v[i][2] = 0.0; + } + } +} + +/* ---------------------------------------------------------------------- + add back in velocity bias to atom I removed by remove_bias() + assume remove_bias() was previously called +------------------------------------------------------------------------- */ + +void ComputeTempPartial::restore_bias(double *v) +{ + v[0] += vbias[0]; + v[1] += vbias[1]; + v[2] += vbias[2]; + if (tbias) tbias->restore_bias(v); +} + +/* ---------------------------------------------------------------------- + add back in velocity bias to all atoms removed by remove_bias_all() + assume remove_bias_all() was previously called +------------------------------------------------------------------------- */ + +void ComputeTempPartial::restore_bias_all() +{ + double **v = atom->v; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + v[i][0] += vbiasall[i][0]; + v[i][1] += vbiasall[i][1]; + v[i][2] += vbiasall[i][2]; + } + + if (tbias) tbias->restore_bias_all(); +} + +/* ---------------------------------------------------------------------- */ + +double ComputeTempPartial::memory_usage() +{ + double bytes = maxbias * sizeof(double); + return bytes; +} diff --git a/src/compute_temp_partial.h b/src/compute_temp_partial.h index b1db6c08cf..a327ce6391 100644 --- a/src/compute_temp_partial.h +++ b/src/compute_temp_partial.h @@ -25,13 +25,23 @@ class ComputeTempPartial : public Compute { void init(); double compute_scalar(); void compute_vector(); + void remove_bias(int, double *); + void remove_bias_all(); + void restore_bias(double *); + void restore_bias_all(); + double memory_usage(); private: int xflag,yflag,zflag; int fix_dof; double tfactor; + double vbias[3]; // stored velocity bias for one atom + double **vbiasall; // stored velocity bias for all atoms + int maxbias; // size of vbiasall array + Compute *tbias; // ptr to additional bias compute + void recount(); }; diff --git a/src/compute_temp_ramp.cpp b/src/compute_temp_ramp.cpp index e88b446a73..91f6409816 100644 --- a/src/compute_temp_ramp.cpp +++ b/src/compute_temp_ramp.cpp @@ -22,6 +22,7 @@ #include "fix.h" #include "domain.h" #include "lattice.h" +#include "memory.h" #include "error.h" using namespace LAMMPS_NS; @@ -108,8 +109,9 @@ ComputeTempRamp::ComputeTempRamp(LAMMPS *lmp, int narg, char **arg) : extvector = 1; tempflag = 1; tempbias = 1; - vbias[0] = vbias[1] = vbias[2] = 0.0; + maxbias = 0; + vbiasall = NULL; vector = new double[6]; } @@ -117,6 +119,7 @@ ComputeTempRamp::ComputeTempRamp(LAMMPS *lmp, int narg, char **arg) : ComputeTempRamp::~ComputeTempRamp() { + memory->destroy_2d_double_array(vbiasall); delete [] vector; } @@ -128,6 +131,12 @@ void ComputeTempRamp::init() for (int i = 0; i < modify->nfix; i++) fix_dof += modify->fix[i]->dof(igroup); recount(); + + if (id_bias) { + int i = modify->find_compute(id_bias); + if (i < 0) error->all("Could not find compute ID for temperature bias"); + tbias = modify->compute[i]; + } } /* ---------------------------------------------------------------------- */ @@ -149,6 +158,12 @@ double ComputeTempRamp::compute_scalar() invoked |= INVOKED_SCALAR; + if (tbias) { + if (!(tbias->invoked & INVOKED_SCALAR)) + double tmp = tbias->compute_scalar(); + tbias->remove_bias_all(); + } + double **x = atom->x; double **v = atom->v; double *mass = atom->mass; @@ -176,6 +191,8 @@ double ComputeTempRamp::compute_scalar() vthermal[2]*vthermal[2]) * rmass[i]; } + if (tbias) tbias->restore_bias_all(); + MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic) recount(); scalar *= tfactor; @@ -191,6 +208,11 @@ void ComputeTempRamp::compute_vector() invoked |= INVOKED_VECTOR; + if (tbias) { + if (!(tbias->invoked & INVOKED_VECTOR)) tbias->compute_vector(); + tbias->remove_bias_all(); + } + double **x = atom->x; double **v = atom->v; double *mass = atom->mass; @@ -223,17 +245,90 @@ void ComputeTempRamp::compute_vector() t[5] += massone * vthermal[1]*vthermal[2]; } + if (tbias) tbias->restore_bias_all(); + MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); for (i = 0; i < 6; i++) vector[i] *= force->mvv2e; } -/* ---------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + remove velocity bias from atom I to leave thermal velocity +------------------------------------------------------------------------- */ void ComputeTempRamp::remove_bias(int i, double *v) { + if (tbias) tbias->remove_bias(i,v); + double fraction = (atom->x[i][coord_dim] - coord_lo) / (coord_hi - coord_lo); fraction = MAX(fraction,0.0); fraction = MIN(fraction,1.0); vbias[v_dim] = v_lo + fraction*(v_hi - v_lo); v[v_dim] -= vbias[v_dim]; } + +/* ---------------------------------------------------------------------- + remove velocity bias from all atoms to leave thermal velocity +------------------------------------------------------------------------- */ + +void ComputeTempRamp::remove_bias_all() +{ + if (tbias) tbias->remove_bias_all(); + + double **v = atom->v; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + if (nlocal > maxbias) { + memory->destroy_2d_double_array(vbiasall); + maxbias = atom->nmax; + vbiasall = memory->create_2d_double_array(maxbias,3, + "compute/temp:vbiasall"); + } + + double fraction; + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + fraction = (atom->x[i][coord_dim] - coord_lo) / (coord_hi - coord_lo); + fraction = MAX(fraction,0.0); + fraction = MIN(fraction,1.0); + vbiasall[i][v_dim] = v_lo + fraction*(v_hi - v_lo); + v[i][v_dim] -= vbiasall[i][v_dim]; + } +} + +/* ---------------------------------------------------------------------- + add back in velocity bias to atom I removed by remove_bias() + assume remove_bias() was previously called +------------------------------------------------------------------------- */ + +void ComputeTempRamp::restore_bias(double *v) +{ + v[v_dim] += vbias[v_dim]; + if (tbias) tbias->restore_bias(v); +} + +/* ---------------------------------------------------------------------- + add back in velocity bias to all atoms removed by remove_bias_all() + assume remove_bias_all() was previously called +------------------------------------------------------------------------- */ + +void ComputeTempRamp::restore_bias_all() +{ + double **v = atom->v; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + v[i][v_dim] += vbiasall[i][v_dim]; + + if (tbias) tbias->restore_bias_all(); +} + +/* ---------------------------------------------------------------------- */ + +double ComputeTempRamp::memory_usage() +{ + double bytes = maxbias * sizeof(double); + return bytes; +} diff --git a/src/compute_temp_ramp.h b/src/compute_temp_ramp.h index 31b624cbc4..44753eeba8 100644 --- a/src/compute_temp_ramp.h +++ b/src/compute_temp_ramp.h @@ -25,7 +25,12 @@ class ComputeTempRamp : public Compute { void init(); double compute_scalar(); void compute_vector(); + void remove_bias(int, double *); + void remove_bias_all(); + void restore_bias(double *); + void restore_bias_all(); + double memory_usage(); private: int coord_dim; @@ -35,6 +40,11 @@ class ComputeTempRamp : public Compute { int scaleflag,fix_dof; double tfactor,xscale,yscale,zscale; + double vbias[3]; // stored velocity bias for one atom + double **vbiasall; // stored velocity bias for all atoms + int maxbias; // size of vbiasall array + Compute *tbias; // ptr to additional bias compute + void recount(); }; diff --git a/src/compute_temp_region.cpp b/src/compute_temp_region.cpp index 3d40dcde6e..efab6219de 100644 --- a/src/compute_temp_region.cpp +++ b/src/compute_temp_region.cpp @@ -16,8 +16,10 @@ #include "compute_temp_region.h" #include "atom.h" #include "force.h" +#include "modify.h" #include "domain.h" #include "region.h" +#include "memory.h" #include "error.h" using namespace LAMMPS_NS; @@ -42,6 +44,8 @@ ComputeTempRegion::ComputeTempRegion(LAMMPS *lmp, int narg, char **arg) : tempflag = 1; tempbias = 1; + maxbias = 0; + vbiasall = NULL; vector = new double[6]; } @@ -49,6 +53,7 @@ ComputeTempRegion::ComputeTempRegion(LAMMPS *lmp, int narg, char **arg) : ComputeTempRegion::~ComputeTempRegion() { + memory->destroy_2d_double_array(vbiasall); delete [] vector; } @@ -57,6 +62,12 @@ ComputeTempRegion::~ComputeTempRegion() void ComputeTempRegion::init() { dof = 0; + + if (id_bias) { + int i = modify->find_compute(id_bias); + if (i < 0) error->all("Could not find compute ID for temperature bias"); + tbias = modify->compute[i]; + } } /* ---------------------------------------------------------------------- */ @@ -65,6 +76,12 @@ double ComputeTempRegion::compute_scalar() { invoked |= INVOKED_SCALAR; + if (tbias) { + if (!(tbias->invoked & INVOKED_SCALAR)) + double tmp = tbias->compute_scalar(); + tbias->remove_bias_all(); + } + double **x = atom->x; double **v = atom->v; double *mass = atom->mass; @@ -93,6 +110,8 @@ double ComputeTempRegion::compute_scalar() } } + if (tbias) tbias->restore_bias_all(); + double tarray[2],tarray_all[2]; tarray[0] = count; tarray[1] = t; @@ -111,6 +130,11 @@ void ComputeTempRegion::compute_vector() invoked |= INVOKED_VECTOR; + if (tbias) { + if (!(tbias->invoked & INVOKED_VECTOR)) tbias->compute_vector(); + tbias->remove_bias_all(); + } + double **x = atom->x; double **v = atom->v; double *mass = atom->mass; @@ -135,14 +159,20 @@ void ComputeTempRegion::compute_vector() t[5] += massone * v[i][1]*v[i][2]; } + if (tbias) tbias->restore_bias_all(); + MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); for (i = 0; i < 6; i++) vector[i] *= force->mvv2e; } -/* ---------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + remove velocity bias from atom I to leave thermal velocity +------------------------------------------------------------------------- */ void ComputeTempRegion::remove_bias(int i, double *v) { + if (tbias) tbias->remove_bias(i,v); + double *x = atom->x[i]; if (atom->mask[i] & groupbit && domain->regions[iregion]->match(x[0],x[1],x[2])) @@ -154,3 +184,79 @@ void ComputeTempRegion::remove_bias(int i, double *v) v[0] = v[1] = v[2] = 0.0; } } + +/* ---------------------------------------------------------------------- + remove velocity bias from all atoms to leave thermal velocity +------------------------------------------------------------------------- */ + +void ComputeTempRegion::remove_bias_all() +{ + if (tbias) tbias->remove_bias_all(); + + double **x = atom->x; + double **v = atom->v; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + if (nlocal > maxbias) { + memory->destroy_2d_double_array(vbiasall); + maxbias = atom->nmax; + vbiasall = memory->create_2d_double_array(maxbias,3, + "compute/temp:vbiasall"); + } + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + if (atom->mask[i] & groupbit && + domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2])) + vbiasall[i][0] = vbiasall[i][1] = vbiasall[i][2] = 0.0; + else { + vbiasall[i][0] = v[i][0]; + vbiasall[i][1] = v[i][1]; + vbiasall[i][2] = v[i][2]; + v[i][0] = v[i][1] = v[i][2] = 0.0; + } + } +} + +/* ---------------------------------------------------------------------- + add back in velocity bias to atom I removed by remove_bias() + assume remove_bias() was previously called +------------------------------------------------------------------------- */ + +void ComputeTempRegion::restore_bias(double *v) +{ + v[0] += vbias[0]; + v[1] += vbias[1]; + v[2] += vbias[2]; + if (tbias) tbias->restore_bias(v); +} + +/* ---------------------------------------------------------------------- + add back in velocity bias to all atoms removed by remove_bias_all() + assume remove_bias_all() was previously called +------------------------------------------------------------------------- */ + +void ComputeTempRegion::restore_bias_all() +{ + double **v = atom->v; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + v[i][0] += vbiasall[i][0]; + v[i][1] += vbiasall[i][1]; + v[i][2] += vbiasall[i][2]; + } + + if (tbias) tbias->restore_bias_all(); +} + +/* ---------------------------------------------------------------------- */ + +double ComputeTempRegion::memory_usage() +{ + double bytes = maxbias * sizeof(double); + return bytes; +} diff --git a/src/compute_temp_region.h b/src/compute_temp_region.h index c8944e2aed..5d434c077b 100644 --- a/src/compute_temp_region.h +++ b/src/compute_temp_region.h @@ -25,10 +25,20 @@ class ComputeTempRegion : public Compute { void init(); double compute_scalar(); void compute_vector(); + void remove_bias(int, double *); + void remove_bias_all(); + void restore_bias(double *); + void restore_bias_all(); + double memory_usage(); private: int iregion; + + double vbias[3]; // stored velocity bias for one atom + double **vbiasall; // stored velocity bias for all atoms + int maxbias; // size of vbiasall array + Compute *tbias; // ptr to additional bias compute }; } diff --git a/src/compute_temp_sphere.cpp b/src/compute_temp_sphere.cpp index 24fb118ca2..4145db51c7 100644 --- a/src/compute_temp_sphere.cpp +++ b/src/compute_temp_sphere.cpp @@ -62,6 +62,13 @@ void ComputeTempSphere::init() fix_dof += modify->fix[i]->dof(igroup); recount(); + if (id_bias) { + tempbias = 1; + int i = modify->find_compute(id_bias); + if (i < 0) error->all("Could not find compute ID for temperature bias"); + tbias = modify->compute[i]; + } else tempbias = 0; + if (atom->mass) { double *mass = atom->mass; double **shape = atom->shape; @@ -92,6 +99,12 @@ double ComputeTempSphere::compute_scalar() { invoked |= INVOKED_SCALAR; + if (tbias) { + if (!(tbias->invoked & INVOKED_SCALAR)) + double tmp = tbias->compute_scalar(); + tbias->remove_bias_all(); + } + double **v = atom->v; double **omega = atom->omega; double *mass = atom->mass; @@ -120,6 +133,8 @@ double ComputeTempSphere::compute_scalar() } } + if (tbias) tbias->restore_bias_all(); + MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic) recount(); scalar *= tfactor; @@ -134,6 +149,11 @@ void ComputeTempSphere::compute_vector() invoked |= INVOKED_VECTOR; + if (tbias) { + if (!(tbias->invoked & INVOKED_VECTOR)) tbias->compute_vector(); + tbias->remove_bias_all(); + } + double **v = atom->v; double **omega = atom->omega; double *mass = atom->mass; @@ -186,6 +206,47 @@ void ComputeTempSphere::compute_vector() } } + if (tbias) tbias->restore_bias_all(); + MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); for (i = 0; i < 6; i++) vector[i] *= force->mvv2e; } + +/* ---------------------------------------------------------------------- + remove velocity bias from atom I to leave thermal velocity +------------------------------------------------------------------------- */ + +void ComputeTempSphere::remove_bias(int i, double *v) +{ + if (tbias) tbias->remove_bias(i,v); +} + +/* ---------------------------------------------------------------------- + remove velocity bias from all atoms to leave thermal velocity +------------------------------------------------------------------------- */ + +void ComputeTempSphere::remove_bias_all() +{ + if (tbias) tbias->remove_bias_all(); +} + +/* ---------------------------------------------------------------------- + add back in velocity bias to atom I removed by remove_bias() + assume remove_bias() was previously called +------------------------------------------------------------------------- */ + +void ComputeTempSphere::restore_bias(double *v) +{ + if (tbias) tbias->restore_bias(v); +} + +/* ---------------------------------------------------------------------- + add back in velocity bias to all atoms removed by remove_bias_all() + assume remove_bias_all() was previously called +------------------------------------------------------------------------- */ + +void ComputeTempSphere::restore_bias_all() +{ + if (tbias) tbias->restore_bias_all(); +} + diff --git a/src/compute_temp_sphere.h b/src/compute_temp_sphere.h index 6ae185f3a7..d056ee7b5d 100644 --- a/src/compute_temp_sphere.h +++ b/src/compute_temp_sphere.h @@ -26,11 +26,18 @@ class ComputeTempSphere : public Compute { double compute_scalar(); void compute_vector(); + void remove_bias(int, double *); + void remove_bias_all(); + void restore_bias(double *); + void restore_bias_all(); + private: int fix_dof; double tfactor; double *inertia; + Compute *tbias; // ptr to additional bias compute + void recount(); }; diff --git a/src/style_asphere.h b/src/style_asphere.h index 1a84677e10..70dba98aa4 100644 --- a/src/style_asphere.h +++ b/src/style_asphere.h @@ -25,7 +25,7 @@ AtomStyle(ellipsoid,AtomVecEllipsoid) #endif #ifdef ComputeClass -ComputeStyle(erotate/asphere,ComputeERotateASphere) +ComputeStyle(erotate/asphere,ComputeERotateAsphere) ComputeStyle(temp/asphere,ComputeTempAsphere) #endif @@ -36,9 +36,9 @@ ComputeStyle(temp/asphere,ComputeTempAsphere) #endif #ifdef FixClass -FixStyle(nve/asphere,FixNVEASphere) -FixStyle(nvt/asphere,FixNVTASphere) -FixStyle(npt/asphere,FixNPTASphere) +FixStyle(nve/asphere,FixNVEAsphere) +FixStyle(nvt/asphere,FixNVTAsphere) +FixStyle(npt/asphere,FixNPTAsphere) #endif #ifdef PairInclude