ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

mods to change_box

Browse Source
This commit is contained in:
Steve Plimpton
2023-09-01 12:15:51 -06:00
parent 932a080246
commit 7462439b5d
6 changed files with 157 additions and 37 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
src/create_atoms.cpp
View File

@ -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) {

134
src/create_box.cpp
View File

@ -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<RegPrism *>(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<RegPrism *>(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);

5
src/domain.cpp
View File

@ -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];

30
src/lattice.cpp
View File

@ -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
------------------------------------------------------------------------- */

10
src/lattice.h
View File

@ -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];

13
src/read_data.cpp
View File

@ -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) {