From 7462439b5d818f7d2fbc8ab74b344f9b76e4eeee Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Fri, 1 Sep 2023 12:15:51 -0600 Subject: [PATCH] mods to change_box --- src/create_atoms.cpp | 2 +- src/create_box.cpp | 134 +++++++++++++++++++++++++++++++++---------- src/domain.cpp | 5 ++ src/lattice.cpp | 30 ++++++++++ src/lattice.h | 10 +++- src/read_data.cpp | 13 ++++- 6 files changed, 157 insertions(+), 37 deletions(-) diff --git a/src/create_atoms.cpp b/src/create_atoms.cpp index 32be85e647..aa95c8397f 100644 --- a/src/create_atoms.cpp +++ b/src/create_atoms.cpp @@ -1178,7 +1178,7 @@ void CreateAtoms::add_lattice() } else domain->bbox(domain->sublo_lamda, domain->subhi_lamda, bboxlo, bboxhi); - // narrow down the subbox by the bounding box of the given region, if available. + // narrow down the subbox by the bounding box of the given region, if available // for small regions in large boxes, this can result in a significant speedup if ((style == REGION) && region->bboxflag) { diff --git a/src/create_box.cpp b/src/create_box.cpp index 02aa63acf0..74a8db8bb3 100644 --- a/src/create_box.cpp +++ b/src/create_box.cpp @@ -19,6 +19,7 @@ #include "domain.h" #include "error.h" #include "force.h" +#include "lattice.h" #include "region.h" #include "region_prism.h" #include "update.h" @@ -45,40 +46,112 @@ void CreateBox::command(int narg, char **arg) // region check - auto region = domain->get_region_by_id(arg[1]); - if (!region) error->all(FLERR, "Create_box region {} does not exist", arg[1]); - if (region->bboxflag == 0) error->all(FLERR, "Create_box region does not support a bounding box"); + Region *region = nullptr; + int triclinic_general = 0; + + if (strcmp(arg[1],"NULL") == 0) triclinic_general = 1; + else { + region = domain->get_region_by_id(arg[1]); + if (!region) error->all(FLERR, "Create_box region {} does not exist", arg[1]); + if (region->bboxflag == 0) error->all(FLERR, "Create_box region does not support a bounding box"); + region->init(); + } - region->init(); + // setup simulation box + // 3 options: orthogonal, restricted triclinic, general triclinic - // if region not prism: - // setup orthogonal domain - // set simulation domain from region extent - // if region is prism: - // seutp triclinic domain - // set simulation domain params from prism params + int iarg = 2; + + if (region) { + // region is not prism + // setup orthogonal box + // set simulation domain from region extent + + if (strcmp(region->style, "prism") != 0) { + domain->triclinic = 0; + domain->boxlo[0] = region->extent_xlo; + domain->boxhi[0] = region->extent_xhi; + domain->boxlo[1] = region->extent_ylo; + domain->boxhi[1] = region->extent_yhi; + domain->boxlo[2] = region->extent_zlo; + domain->boxhi[2] = region->extent_zhi; - if (strcmp(region->style, "prism") != 0) { - domain->triclinic = 0; - domain->boxlo[0] = region->extent_xlo; - domain->boxhi[0] = region->extent_xhi; - domain->boxlo[1] = region->extent_ylo; - domain->boxhi[1] = region->extent_yhi; - domain->boxlo[2] = region->extent_zlo; - domain->boxhi[2] = region->extent_zhi; + // region is prism + // seutp restricted triclinic box + // set simulation domain from prism params - } else { - domain->triclinic = 1; - auto prism = dynamic_cast(region); - domain->boxlo[0] = prism->xlo; - domain->boxhi[0] = prism->xhi; - domain->boxlo[1] = prism->ylo; - domain->boxhi[1] = prism->yhi; - domain->boxlo[2] = prism->zlo; - domain->boxhi[2] = prism->zhi; - domain->xy = prism->xy; - domain->xz = prism->xz; - domain->yz = prism->yz; + } else { + domain->triclinic = 1; + auto prism = dynamic_cast(region); + domain->boxlo[0] = prism->xlo; + domain->boxhi[0] = prism->xhi; + domain->boxlo[1] = prism->ylo; + domain->boxhi[1] = prism->yhi; + domain->boxlo[2] = prism->zlo; + domain->boxhi[2] = prism->zhi; + domain->xy = prism->xy; + domain->xz = prism->xz; + domain->yz = prism->yz; + } + + // setup general triclinic box (with no region) + // read next box extent arguments to create ABC edge vectors + origin + // set_general_triclinic() converts + // ABC edge vectors + origin to restricted triclinic + + } else if (triclinic_general) { + if (!domain->lattice->is_custom()) + error->all(FLERR,"Create_box with no region requires custom lattice"); + if (domain->lattice->is_oriented()) + error->all(FLERR,"Create_box with no region requires lattice with default orientation"); + + if (iarg + 6 > narg) utils::missing_cmd_args(FLERR, "create_box general triclinic", error); + + double alo = utils::numeric(FLERR, arg[iarg + 0], false, lmp); + double ahi = utils::numeric(FLERR, arg[iarg + 1], false, lmp); + double blo = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + double bhi = utils::numeric(FLERR, arg[iarg + 3], false, lmp); + double clo = utils::numeric(FLERR, arg[iarg + 4], false, lmp); + double chi = utils::numeric(FLERR, arg[iarg + 5], false, lmp); + iarg += 6; + + // use lattice2box() to generate avec, bvec, cvec and origin + + double avec[3],bvec[3],cvec[3],origin[3]; + double px,py,pz; + + px = alo; py = blo; pz = clo; + domain->lattice->lattice2box(px,py,pz); + origin[0] = px; + origin[1] = py; + origin[2] = pz; + + px = ahi; py = blo; pz = clo; + domain->lattice->lattice2box(px,py,pz); + avec[0] = px - origin[0]; + avec[1] = py - origin[1]; + avec[2] = pz - origin[2]; + + px = alo; py = bhi; pz = clo; + domain->lattice->lattice2box(px,py,pz); + bvec[0] = px - origin[0]; + bvec[1] = py - origin[1]; + bvec[2] = pz - origin[2]; + + px = alo; py = blo; pz = chi; + domain->lattice->lattice2box(px,py,pz); + cvec[0] = px - origin[0]; + cvec[1] = py - origin[1]; + cvec[2] = pz - origin[2]; + + printf("ORIGIN %g %g %g\n",origin[0],origin[1],origin[2]); + printf("AVEC %g %g %g\n",avec[0],avec[1],avec[2]); + printf("BVEC %g %g %g\n",bvec[0],bvec[1],bvec[2]); + printf("CVEC %g %g %g\n",cvec[0],cvec[1],cvec[2]); + + // setup general triclinic box in Domain + + domain->set_general_triclinic(avec,bvec,cvec,origin); } // if molecular, zero out topology info @@ -104,7 +177,6 @@ void CreateBox::command(int narg, char **arg) // process optional args that can overwrite default settings - int iarg = 2; while (iarg < narg) { if (strcmp(arg[iarg], "bond/types") == 0) { if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "create_box bond/type", error); diff --git a/src/domain.cpp b/src/domain.cpp index 6523bfc8c7..9629645aa6 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -551,6 +551,11 @@ void Domain::set_general_triclinic(double *avec_caller, double *bvec_caller, tri_origin[1] = origin_caller[1]; tri_origin[2] = origin_caller[2]; + // error check on cvec for 2d systems + + if (dimension == 2 && (cvec[0] != 0.0 || cvec[1] != 0.0)) + error->all(FLERR,"General triclinic box edge vector C invalid for 2d system"); + // error check for co-planar A,B,C double abcross[3]; diff --git a/src/lattice.cpp b/src/lattice.cpp index edb482cfac..539ec73531 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -262,6 +262,14 @@ Lattice::Lattice(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) scale = pow(nbasis/volume/scale,1.0/dimension); } + // set orientflag + // for general triclinic, create_box and create_atoms require orientflag = 0 + + oriented = 0; + if (orientx[0] != 1 || orientx[1] != 0 || orientx[2] != 0) oriented = 1; + if (orienty[0] != 0 || orienty[1] != 1 || orienty[2] != 0) oriented = 1; + if (orientz[0] != 0 || orientz[1] != 0 || orientz[2] != 1) oriented = 1; + // initialize lattice <-> box transformation matrices setup_transform(); @@ -311,6 +319,28 @@ Lattice::~Lattice() memory->destroy(basis); } +/* ---------------------------------------------------------------------- + return 1 if style = CUSTOM, else 0 + queried by create_box and create_atoms for general triclinic +------------------------------------------------------------------------- */ + +int Lattice::is_custom() +{ + if (style == CUSTOM) return 1; + return 0; +} + +/* ---------------------------------------------------------------------- + return 1 any orient vectors are non-default, else 0 + queried by create_box and create_atoms for general triclinic +------------------------------------------------------------------------- */ + +int Lattice::is_oriented() +{ + if (oriented) return 1; + return 0; +} + /* ---------------------------------------------------------------------- check if 3 orientation vectors are mutually orthogonal ------------------------------------------------------------------------- */ diff --git a/src/lattice.h b/src/lattice.h index 5b98f580b7..d91b5cc834 100644 --- a/src/lattice.h +++ b/src/lattice.h @@ -33,15 +33,19 @@ class Lattice : protected Pointers { ~Lattice() override; void lattice2box(double &, double &, double &); void box2lattice(double &, double &, double &); - void bbox(int, double, double, double, double &, double &, double &, double &, double &, - double &); - + void bbox(int, double, double, double, + double &, double &, double &, double &, double &, double &); + int is_custom(); + int is_oriented(); + private: double scale; double origin[3]; // lattice origin + int oriented; // 1 if non-default orient xyz, else 0 int orientx[3]; // lattice orientation vecs int orienty[3]; // orientx = what lattice dir lies int orientz[3]; // along x dim in box + double primitive[3][3]; // lattice <-> box transform matrices double priminv[3][3]; diff --git a/src/read_data.cpp b/src/read_data.cpp index 8df901062c..bd998ed3f5 100644 --- a/src/read_data.cpp +++ b/src/read_data.cpp @@ -557,9 +557,11 @@ void ReadData::command(int narg, char **arg) // setup simulation box // 3 options: orthogonal, restricted triclinic, general triclinic - // for general triclinic: convect general ABC edge vectors to LAMMPS restricted triclinic if (!triclinic_general) { + + // orthongal box + domain->boxlo[0] = boxlo[0]; domain->boxhi[0] = boxhi[0]; domain->boxlo[1] = boxlo[1]; @@ -567,6 +569,8 @@ void ReadData::command(int narg, char **arg) domain->boxlo[2] = boxlo[2]; domain->boxhi[2] = boxhi[2]; + // restricted triclinic box + if (triclinic) { domain->triclinic = 1; domain->xy = xy; @@ -574,6 +578,10 @@ void ReadData::command(int narg, char **arg) domain->yz = yz; } + // general triclinic box + // set_general_triclinic() converts + // ABC edge vectors + origin to restricted triclinic + } else if (triclinic_general) { domain->set_general_triclinic(avec,bvec,cvec,tri_origin); } @@ -581,7 +589,8 @@ void ReadData::command(int narg, char **arg) // change simulation box to be union of existing box and new box + shift // only done if firstpass and not first data file - // shift not allowed for general triclinic + // for restricted triclinic, new tilt factors not allowed + // for general triclinic, different new box and shift not allowed if (firstpass && addflag != NONE) {