diff --git a/src/COLLOID/pair_lubricate.cpp b/src/COLLOID/pair_lubricate.cpp index 78c7088fa8..1c809f0f7f 100644 --- a/src/COLLOID/pair_lubricate.cpp +++ b/src/COLLOID/pair_lubricate.cpp @@ -86,6 +86,7 @@ void PairLubricate::compute(int eflag, int vflag) int nlocal = atom->nlocal; int newton_pair = force->newton_pair; int omega_flag = atom->omega_flag; + double vxmu2f = force->vxmu2f; inum = list->inum; ilist = list->ilist; @@ -198,16 +199,17 @@ void PairLubricate::compute(int eflag, int vflag) h_sep = r - 2.0*radi; if (flag1) - a_squeeze = (3.0*PI*mu*2.0*radi/2.0) * (2.0*radi/4.0/h_sep); + a_squeeze = vxmu2f * + (3.0*PI*mu*2.0*radi/2.0) * (2.0*radi/4.0/h_sep); if (flag2) - a_shear = (PI*mu*2.*radi/2.0) * + a_shear = vxmu2f * (PI*mu*2.*radi/2.0) * log(2.0*radi/2.0/h_sep)*(2.0*radi+h_sep)*(2.0*radi+h_sep)/4.0; if (flag3) - a_pump = (PI*mu*pow(2.0*radi,4)/8.0) * + a_pump = vxmu2f * (PI*mu*pow(2.0*radi,4)/8.0) * ((3.0/20.0) * log(2.0*radi/2.0/h_sep) + (63.0/250.0) * (h_sep/2.0/radi) * log(2.0*radi/2.0/h_sep)); if (flag4) - a_twist = (PI*mu*pow(2.0*radi,4)/4.0) * + a_twist = vxmu2f * (PI*mu*pow(2.0*radi,4)/4.0) * (h_sep/2.0/radi) * log(2.0/(2.0*h_sep)); if (h_sep >= cut_inner[itype][jtype]) { diff --git a/src/fix_nph.cpp b/src/fix_nph.cpp index d5f489b3ee..98b8fbb2f1 100644 --- a/src/fix_nph.cpp +++ b/src/fix_nph.cpp @@ -37,6 +37,8 @@ using namespace LAMMPS_NS; #define MIN(A,B) ((A) < (B)) ? (A) : (B) #define MAX(A,B) ((A) > (B)) ? (A) : (B) +enum{XYZ,XY,YZ,XZ,ANISO}; + /* ---------------------------------------------------------------------- */ FixNPH::FixNPH(LAMMPS *lmp, int narg, char **arg) : @@ -55,7 +57,7 @@ FixNPH::FixNPH(LAMMPS *lmp, int narg, char **arg) : if (strcmp(arg[3],"xyz") == 0) { if (narg < 7) error->all("Illegal fix nph command"); - press_couple = 0; + press_couple = XYZ; p_start[0] = p_start[1] = p_start[2] = atof(arg[4]); p_stop[0] = p_stop[1] = p_stop[2] = atof(arg[5]); p_period[0] = p_period[1] = p_period[2] = atof(arg[6]); @@ -66,16 +68,16 @@ FixNPH::FixNPH(LAMMPS *lmp, int narg, char **arg) : } } else { - if (strcmp(arg[3],"xy") == 0) press_couple = 1; - else if (strcmp(arg[3],"yz") == 0) press_couple = 2; - else if (strcmp(arg[3],"xz") == 0) press_couple = 3; - else if (strcmp(arg[3],"aniso") == 0) press_couple = 4; + if (strcmp(arg[3],"xy") == 0) press_couple = XY; + else if (strcmp(arg[3],"yz") == 0) press_couple = YZ; + else if (strcmp(arg[3],"xz") == 0) press_couple = XZ; + else if (strcmp(arg[3],"aniso") == 0) press_couple = ANISO; else error->all("Illegal fix nph command"); if (narg < 11) error->all("Illegal fix nph command"); if (domain->dimension == 2 && - (press_couple == 1 || press_couple == 2 || press_couple == 3)) + (press_couple == XY || press_couple == YZ || press_couple == XZ)) error->all("Invalid fix nph command for a 2d simulation"); if (strcmp(arg[4],"NULL") == 0) { @@ -117,7 +119,7 @@ FixNPH::FixNPH(LAMMPS *lmp, int narg, char **arg) : allremap = 1; int iarg; - if (press_couple == 0) iarg = 7; + if (press_couple == XYZ) iarg = 7; else iarg = 11; while (iarg < narg) { @@ -134,7 +136,7 @@ FixNPH::FixNPH(LAMMPS *lmp, int narg, char **arg) : } else error->all("Illegal fix nph command"); } - // error checks + // check for periodicity in controlled dimensions if (p_flag[0] && domain->xperiodic == 0) error->all("Cannot use fix nph on a non-periodic dimension"); @@ -313,7 +315,7 @@ void FixNPH::setup(int vflag) p_target[1] = p_start[1]; p_target[2] = p_start[2]; - if (press_couple == 0) { + if (press_couple == XYZ) { double tmp = temperature->compute_scalar(); tmp = pressure->compute_scalar(); } else { @@ -426,7 +428,7 @@ void FixNPH::final_integrate() // compute new pressure - if (press_couple == 0) { + if (press_couple == XYZ) { double tmp = temperature->compute_scalar(); tmp = pressure->compute_scalar(); } else { @@ -598,21 +600,21 @@ void FixNPH::couple() { double *tensor = pressure->vector; - if (press_couple == 0) + if (press_couple == XYZ) p_current[0] = p_current[1] = p_current[2] = pressure->scalar; - else if (press_couple == 1) { + else if (press_couple == XY) { double ave = 0.5 * (tensor[0] + tensor[1]); p_current[0] = p_current[1] = ave; p_current[2] = tensor[2]; - } else if (press_couple == 2) { + } else if (press_couple == YZ) { double ave = 0.5 * (tensor[1] + tensor[2]); p_current[1] = p_current[2] = ave; p_current[0] = tensor[0]; - } else if (press_couple == 3) { + } else if (press_couple == XZ) { double ave = 0.5 * (tensor[0] + tensor[2]); p_current[0] = p_current[2] = ave; p_current[1] = tensor[1]; - } if (press_couple == 4) { + } else if (press_couple == ANISO) { p_current[0] = tensor[0]; p_current[1] = tensor[1]; p_current[2] = tensor[2]; diff --git a/src/fix_npt.cpp b/src/fix_npt.cpp index c3c945ba44..04d3efb30d 100644 --- a/src/fix_npt.cpp +++ b/src/fix_npt.cpp @@ -37,6 +37,7 @@ using namespace LAMMPS_NS; #define MAX(A,B) ((A) > (B)) ? (A) : (B) enum{NOBIAS,BIAS}; +enum{XYZ,XY,YZ,XZ,ANISO}; /* ---------------------------------------------------------------------- */ @@ -63,7 +64,7 @@ FixNPT::FixNPT(LAMMPS *lmp, int narg, char **arg) : if (strcmp(arg[6],"xyz") == 0) { if (narg < 10) error->all("Illegal fix npt command"); - press_couple = 0; + press_couple = XYZ; p_start[0] = p_start[1] = p_start[2] = atof(arg[7]); p_stop[0] = p_stop[1] = p_stop[2] = atof(arg[8]); p_period[0] = p_period[1] = p_period[2] = atof(arg[9]); @@ -74,16 +75,16 @@ FixNPT::FixNPT(LAMMPS *lmp, int narg, char **arg) : } } else { - if (strcmp(arg[6],"xy") == 0) press_couple = 1; - else if (strcmp(arg[6],"yz") == 0) press_couple = 2; - else if (strcmp(arg[6],"xz") == 0) press_couple = 3; - else if (strcmp(arg[6],"aniso") == 0) press_couple = 4; + if (strcmp(arg[6],"xy") == 0) press_couple = XY; + else if (strcmp(arg[6],"yz") == 0) press_couple = YZ; + else if (strcmp(arg[6],"xz") == 0) press_couple = XZ; + else if (strcmp(arg[6],"aniso") == 0) press_couple = ANISO; else error->all("Illegal fix npt command"); if (narg < 14) error->all("Illegal fix npt command"); if (domain->dimension == 2 && - (press_couple == 1 || press_couple == 2 || press_couple == 3)) + (press_couple == XY || press_couple == YZ || press_couple == XZ)) error->all("Invalid fix npt command for a 2d simulation"); if (strcmp(arg[7],"NULL") == 0) { @@ -125,7 +126,7 @@ FixNPT::FixNPT(LAMMPS *lmp, int narg, char **arg) : allremap = 1; int iarg; - if (press_couple == 0) iarg = 10; + if (press_couple == XYZ) iarg = 10; else iarg = 14; while (iarg < narg) { @@ -142,7 +143,7 @@ FixNPT::FixNPT(LAMMPS *lmp, int narg, char **arg) : } else error->all("Illegal fix npt command"); } - // error checks + // check for periodicity in controlled dimensions if (p_flag[0] && domain->xperiodic == 0) error->all("Cannot use fix npt on a non-periodic dimension"); @@ -322,7 +323,7 @@ void FixNPT::setup(int vflag) p_target[2] = p_start[2]; t_current = temperature->compute_scalar(); - if (press_couple == 0) { + if (press_couple == XYZ) { double tmp = pressure->compute_scalar(); } else { temperature->compute_vector(); @@ -473,7 +474,7 @@ void FixNPT::final_integrate() // compute new T,P t_current = temperature->compute_scalar(); - if (press_couple == 0) { + if (press_couple == XYZ) { double tmp = pressure->compute_scalar(); } else { temperature->compute_vector(); @@ -687,21 +688,21 @@ void FixNPT::couple() { double *tensor = pressure->vector; - if (press_couple == 0) + if (press_couple == XYZ) p_current[0] = p_current[1] = p_current[2] = pressure->scalar; - else if (press_couple == 1) { + else if (press_couple == XY) { double ave = 0.5 * (tensor[0] + tensor[1]); p_current[0] = p_current[1] = ave; p_current[2] = tensor[2]; - } else if (press_couple == 2) { + } else if (press_couple == YZ) { double ave = 0.5 * (tensor[1] + tensor[2]); p_current[1] = p_current[2] = ave; p_current[0] = tensor[0]; - } else if (press_couple == 3) { + } else if (press_couple == XZ) { double ave = 0.5 * (tensor[0] + tensor[2]); p_current[0] = p_current[2] = ave; p_current[1] = tensor[1]; - } if (press_couple == 4) { + } else if (press_couple == ANISO) { p_current[0] = tensor[0]; p_current[1] = tensor[1]; p_current[2] = tensor[2]; diff --git a/src/force.h b/src/force.h index 0b2b900418..6bbb566488 100644 --- a/src/force.h +++ b/src/force.h @@ -26,6 +26,7 @@ class Force : protected Pointers { double nktv2p; // conversion of NkT/V to pressure double qqr2e; // conversion of q^2/r to energy double qe2f; // conversion of qE to force + double vxmu2f; // conversion of vx mu to force double dielectric; // dielectric constant double qqrd2e; // q^2/r to energy w/ dielectric constant diff --git a/src/min_cg.cpp b/src/min_cg.cpp index 853791ef39..7723729068 100644 --- a/src/min_cg.cpp +++ b/src/min_cg.cpp @@ -186,6 +186,7 @@ void MinCG::run() eng_force(&ntmp,&xtmp,&htmp,&etmp); output->write(update->ntimestep); } + timer->barrier_stop(TIME_LOOP); // delete fix at end of run, so its atom arrays won't persist diff --git a/src/update.cpp b/src/update.cpp index 06816edaf3..9cbdcc3b94 100644 --- a/src/update.cpp +++ b/src/update.cpp @@ -102,6 +102,7 @@ void Update::set_units(const char *style) force->nktv2p = 1.0; force->qqr2e = 1.0; force->qe2f = 1.0; + force->vxmu2f = 1.0; dt = 0.005; neighbor->skin = 0.3; @@ -112,6 +113,7 @@ void Update::set_units(const char *style) force->nktv2p = 68568.415; force->qqr2e = 332.06371; force->qe2f = 23.060549; + force->vxmu2f = 1.4393264316e4; dt = 1.0; neighbor->skin = 2.0; @@ -122,6 +124,7 @@ void Update::set_units(const char *style) force->nktv2p = 1.6021765e6; force->qqr2e = 14.399645; force->qe2f = 1.0; + force->vxmu2f = 0.6241509647; dt = 0.001; neighbor->skin = 2.0; } else error->all("Illegal units command");