From be1e06897bd6a346efae045c083ef257ad536c57 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Tue, 7 Aug 2012 16:12:48 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@8551 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/domain.cpp | 43 ++++++++++++++++++++++++++++++++++--------- src/domain.h | 2 ++ src/fix_deform.cpp | 30 +++++++++++++++++++++--------- src/fix_deform.h | 2 +- src/fix_nh.cpp | 13 ++++++++++--- src/fix_nh.h | 3 ++- src/input.cpp | 10 ++++++++++ src/input.h | 1 + 8 files changed, 81 insertions(+), 23 deletions(-) diff --git a/src/domain.cpp b/src/domain.cpp index ebf34622a9..69da3f2474 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -65,6 +65,8 @@ Domain::Domain(LAMMPS *lmp) : Pointers(lmp) boundary[2][0] = boundary[2][1] = 0; triclinic = 0; + tiltsmall = 1; + boxlo[0] = boxlo[1] = boxlo[2] = -0.5; boxhi[0] = boxhi[1] = boxhi[2] = 0.5; xy = xz = yz = 0.0; @@ -135,17 +137,20 @@ void Domain::set_initial_box() if (boxlo[0] >= boxhi[0] || boxlo[1] >= boxhi[1] || boxlo[2] >= boxhi[2]) error->one(FLERR,"Box bounds are invalid"); - // error check on triclinic tilt factors + if (domain->dimension == 2 && (xz != 0.0 || yz != 0.0)) + error->all(FLERR,"Cannot skew triclinic box in z for 2d simulation"); + + // error check or warning on triclinic tilt factors if (triclinic) { - if (domain->dimension == 2 && (xz != 0.0 || yz != 0.0)) - error->all(FLERR,"Cannot skew triclinic box in z for 2d simulation"); - if (fabs(xy/(boxhi[0]-boxlo[0])) > 0.5) - error->all(FLERR,"Triclinic box skew is too large"); - if (fabs(xz/(boxhi[0]-boxlo[0])) > 0.5) - error->all(FLERR,"Triclinic box skew is too large"); - if (fabs(yz/(boxhi[1]-boxlo[1])) > 0.5) - error->all(FLERR,"Triclinic box skew is too large"); + if ((fabs(xy/(boxhi[0]-boxlo[0])) > 0.5) || + (fabs(xz/(boxhi[0]-boxlo[0])) > 0.5) || + (fabs(yz/(boxhi[1]-boxlo[1])) > 0.5)) { + if (tiltsmall) + error->all(FLERR,"Triclinic box skew is too large"); + else if (comm->me == 0) + error->warning(FLERR,"Triclinic box skew is large"); + } } // set small based on box size and SMALL @@ -1277,6 +1282,26 @@ void Domain::set_boundary(int narg, char **arg, int flag) } } +/* ---------------------------------------------------------------------- + set domain attributes +------------------------------------------------------------------------- */ + +void Domain::set_box(int narg, char **arg) +{ + if (narg < 1) error->all(FLERR,"Illegal box command"); + + int iarg = 0; + while (iarg < narg) { + if (strcmp(arg[iarg],"tilt") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal box command"); + if (strcmp(arg[iarg+1],"small") == 0) tiltsmall = 1; + else if (strcmp(arg[iarg+1],"large") == 0) tiltsmall = 0; + else error->all(FLERR,"Illegal box command"); + iarg += 2; + } else error->all(FLERR,"Illegal box command"); + } +} + /* ---------------------------------------------------------------------- print box info, orthogonal or triclinic ------------------------------------------------------------------------- */ diff --git a/src/domain.h b/src/domain.h index de645fd92d..ccfbd3a39d 100644 --- a/src/domain.h +++ b/src/domain.h @@ -36,6 +36,7 @@ class Domain : protected Pointers { // 3 = shrink-wrap non-per w/ min int triclinic; // 0 = orthog box, 1 = triclinic + int tiltsmall; // 1 if limit tilt, else 0 // orthogonal box double xprd,yprd,zprd; // global box dimensions @@ -111,6 +112,7 @@ class Domain : protected Pointers { void delete_region(int, char **); int find_region(char *); void set_boundary(int, char **, int); + void set_box(int, char **); void print_box(const char *); void boundary_string(char *); diff --git a/src/fix_deform.cpp b/src/fix_deform.cpp index d972d0057e..bc7b5ddfae 100644 --- a/src/fix_deform.cpp +++ b/src/fix_deform.cpp @@ -354,7 +354,8 @@ FixDeform::FixDeform(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) // reneighboring only forced if flips can occur due to shape changes - if (set[3].style || set[4].style || set[5].style) force_reneighbor = 1; + if (flipflag && (set[3].style || set[4].style || set[5].style)) + force_reneighbor = 1; next_reneighbor = -1; nrigid = 0; @@ -565,23 +566,27 @@ void FixDeform::init() // an integer combination of the edge vectors of the unflipped shape matrix // VARIABLE for yz is error, since no way to calculate if box flip occurs // WIGGLE lo/hi flip test is on min/max oscillation limit, not tilt_stop + // only trigger actual errors if flipflag is set if (set[3].style && set[5].style) { int flag = 0; double lo,hi; - if (set[3].style == VARIABLE) + if (flipflag && set[3].style == VARIABLE) error->all(FLERR,"Fix deform cannot use yz variable with xy"); if (set[3].style == WIGGLE) { lo = set[3].tilt_min; hi = set[3].tilt_max; } else lo = hi = set[3].tilt_stop; - if (lo/(set[1].hi_start-set[1].lo_start) < -0.5 || - hi/(set[1].hi_start-set[1].lo_start) > 0.5) flag = 1; - if (set[1].style) { - if (lo/(set[1].hi_stop-set[1].lo_stop) < -0.5 || - hi/(set[1].hi_stop-set[1].lo_stop) > 0.5) flag = 1; + if (flipflag) { + if (lo/(set[1].hi_start-set[1].lo_start) < -0.5 || + hi/(set[1].hi_start-set[1].lo_start) > 0.5) flag = 1; + if (set[1].style) { + if (lo/(set[1].hi_stop-set[1].lo_stop) < -0.5 || + hi/(set[1].hi_stop-set[1].lo_stop) > 0.5) flag = 1; + } + if (flag) + error->all(FLERR,"Fix deform is changing yz too much with xy"); } - if (flag) error->all(FLERR,"Fix deform is changing yz too much with xy"); } // set domain->h_rate values for use by domain and other fixes/computes @@ -838,7 +843,7 @@ void FixDeform::end_of_step() // check yz first since it may change xz, then xz check comes after // flip is performed on next timestep, before reneighboring in pre-exchange() - if (triclinic) { + if (triclinic && flipflag) { double xprd = set[0].hi_target - set[0].lo_target; double yprd = set[1].hi_target - set[1].lo_target; double xprdinv = 1.0 / xprd; @@ -960,6 +965,7 @@ void FixDeform::options(int narg, char **arg) remapflag = X_REMAP; scaleflag = 1; + flipflag = 1; int iarg = 0; while (iarg < narg) { @@ -976,6 +982,12 @@ void FixDeform::options(int narg, char **arg) else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1; else error->all(FLERR,"Illegal fix deform command"); iarg += 2; + } else if (strcmp(arg[iarg],"flip") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix deform command"); + if (strcmp(arg[iarg+1],"yes") == 0) flipflag = 1; + else if (strcmp(arg[iarg+1],"no") == 0) flipflag = 0; + else error->all(FLERR,"Illegal fix deform command"); + iarg += 2; } else error->all(FLERR,"Illegal fix deform command"); } } diff --git a/src/fix_deform.h b/src/fix_deform.h index d4ee95be97..4b23af6df5 100644 --- a/src/fix_deform.h +++ b/src/fix_deform.h @@ -37,7 +37,7 @@ class FixDeform : public Fix { void end_of_step(); private: - int triclinic,scaleflag; + int triclinic,scaleflag,flipflag; int flip,flipxy,flipxz,flipyz; double *h_rate,*h_ratelo; int varflag; // 1 if VARIABLE option is used, 0 if not diff --git a/src/fix_nh.cpp b/src/fix_nh.cpp index fe1e0fee19..3d2ca8aa29 100644 --- a/src/fix_nh.cpp +++ b/src/fix_nh.cpp @@ -75,6 +75,7 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) eta_mass_flag = 1; omega_mass_flag = 0; etap_mass_flag = 0; + flipflag = 1; // turn on tilt factor scaling, whenever applicable @@ -308,6 +309,12 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) else if (strcmp(arg[iarg+1],"no") == 0) scaleyz = 0; else error->all(FLERR,"Illegal fix nvt/npt/nph command"); iarg += 2; + } else if (strcmp(arg[iarg],"flip") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); + if (strcmp(arg[iarg+1],"yes") == 0) flipflag = 1; + else if (strcmp(arg[iarg+1],"no") == 0) flipflag = 0; + else error->all(FLERR,"Illegal fix nvt/npt/nph command"); + iarg += 2; } else if (strcmp(arg[iarg],"fixedpoint") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); fixedpoint[0] = atof(arg[iarg+1]); @@ -437,8 +444,8 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) // reneighboring only forced if flips can occur due to shape changes - if (p_flag[3] || p_flag[4] || p_flag[5]) force_reneighbor = 1; - if (domain->yz != 0.0 || domain->xz != 0.0 || domain->xy != 0.0) + if (flipflag && (p_flag[3] || p_flag[4] || p_flag[5])) force_reneighbor = 1; + if (flipflag && (domain->yz != 0.0 || domain->xz != 0.0 || domain->xy != 0.0)) force_reneighbor = 1; // convert input periods to frequencies @@ -2199,7 +2206,7 @@ void FixNH::pre_exchange() double xprd = domain->xprd; double yprd = domain->yprd; - // flip is triggered when tilt exceeds 0.5 by an amount DELTAFLIP + // flip is only triggered when tilt exceeds 0.5 by DELTAFLIP // this avoids immediate re-flipping due to tilt oscillations double xtiltmax = (0.5+DELTAFLIP)*xprd; diff --git a/src/fix_nh.h b/src/fix_nh.h index 9ab7b51858..fb141bbb7a 100644 --- a/src/fix_nh.h +++ b/src/fix_nh.h @@ -112,8 +112,9 @@ class FixNH : public Fix { int scaleyz; // 1 if yz scaled with lz int scalexz; // 1 if xz scaled with lz int scalexy; // 1 if xy scaled with ly + int flipflag; // 1 if box flips are invoked as needed - double fixedpoint[3]; // Location of dilation fixed-point + double fixedpoint[3]; // location of dilation fixed-point void couple(); void remap(); diff --git a/src/input.cpp b/src/input.cpp index 399499a5c2..1d4a4af560 100644 --- a/src/input.cpp +++ b/src/input.cpp @@ -450,6 +450,7 @@ int Input::execute_command() else if (!strcmp(command,"bond_coeff")) bond_coeff(); else if (!strcmp(command,"bond_style")) bond_style(); else if (!strcmp(command,"boundary")) boundary(); + else if (!strcmp(command,"box")) box(); else if (!strcmp(command,"communicate")) communicate(); else if (!strcmp(command,"compute")) compute(); else if (!strcmp(command,"compute_modify")) compute_modify(); @@ -954,6 +955,15 @@ void Input::boundary() /* ---------------------------------------------------------------------- */ +void Input::box() +{ + if (domain->box_exist) + error->all(FLERR,"Box command after simulation box is defined"); + domain->set_box(narg,arg); +} + +/* ---------------------------------------------------------------------- */ + void Input::communicate() { comm->set(narg,arg); diff --git a/src/input.h b/src/input.h index ca1f0a56b3..fb0fce2d15 100644 --- a/src/input.h +++ b/src/input.h @@ -71,6 +71,7 @@ class Input : protected Pointers { void bond_coeff(); void bond_style(); void boundary(); + void box(); void communicate(); void compute(); void compute_modify();