diff --git a/src/create_atoms.cpp b/src/create_atoms.cpp index aa95c8397f..9a34544627 100644 --- a/src/create_atoms.cpp +++ b/src/create_atoms.cpp @@ -450,8 +450,8 @@ void CreateAtoms::command(int narg, char **arg) atom->nghost = 0; atom->avec->clear_bonus(); - // add atoms/molecules in one of 3 ways - + // add atoms/molecules with appropriate add() method + bigint natoms_previous = atom->natoms; int nlocal_previous = atom->nlocal; @@ -1155,12 +1155,34 @@ void CreateAtoms::add_mesh(const char *filename) } } +/* ---------------------------------------------------------------------- + add many atoms to general triclinic box by looping over lattice +------------------------------------------------------------------------- */ + +void CreateAtoms::add_lattice_triclinic_general() +{ + + +} + /* ---------------------------------------------------------------------- add many atoms by looping over lattice ------------------------------------------------------------------------- */ void CreateAtoms::add_lattice() { + // add atoms on general triclinic lattice if Domain has setting for it + // verify lattice is valid for general triclinic + + int triclinic_general = domain->triclinic_general; + + if (triclinic_general) { + if (!domain->lattice->is_custom()) + error->all(FLERR,"Create_atoms for general triclinic requires custom lattice"); + if (domain->lattice->is_oriented()) + error->all(FLERR,"Create_atoms for general triclinic requires lattice with default orientation"); + } + // convert 8 corners of my subdomain from box coords to lattice coords // for orthogonal, use corner pts of my subbox // for triclinic, use bounding box of my subbox @@ -1202,17 +1224,56 @@ void CreateAtoms::add_lattice() xmin = ymin = zmin = BIG; xmax = ymax = zmax = -BIG; - // convert to lattice coordinates and set bounding box + // convert 8 corner points of bounding box to lattice coordinates + // compute new bounding box (xyz min/max) in lattice coords + // for orthogonal or restricted triclinic, use 8 corner points of bbox lo/hi + + if (!triclinic_general) { + domain->lattice->bbox(1, bboxlo[0], bboxlo[1], bboxlo[2], xmin, ymin, zmin, xmax, ymax, zmax); + domain->lattice->bbox(1, bboxhi[0], bboxlo[1], bboxlo[2], xmin, ymin, zmin, xmax, ymax, zmax); + domain->lattice->bbox(1, bboxlo[0], bboxhi[1], bboxlo[2], xmin, ymin, zmin, xmax, ymax, zmax); + domain->lattice->bbox(1, bboxhi[0], bboxhi[1], bboxlo[2], xmin, ymin, zmin, xmax, ymax, zmax); + domain->lattice->bbox(1, bboxlo[0], bboxlo[1], bboxhi[2], xmin, ymin, zmin, xmax, ymax, zmax); + domain->lattice->bbox(1, bboxhi[0], bboxlo[1], bboxhi[2], xmin, ymin, zmin, xmax, ymax, zmax); + domain->lattice->bbox(1, bboxlo[0], bboxhi[1], bboxhi[2], xmin, ymin, zmin, xmax, ymax, zmax); + domain->lattice->bbox(1, bboxhi[0], bboxhi[1], bboxhi[2], xmin, ymin, zmin, xmax, ymax, zmax); - domain->lattice->bbox(1, bboxlo[0], bboxlo[1], bboxlo[2], xmin, ymin, zmin, xmax, ymax, zmax); - domain->lattice->bbox(1, bboxhi[0], bboxlo[1], bboxlo[2], xmin, ymin, zmin, xmax, ymax, zmax); - domain->lattice->bbox(1, bboxlo[0], bboxhi[1], bboxlo[2], xmin, ymin, zmin, xmax, ymax, zmax); - domain->lattice->bbox(1, bboxhi[0], bboxhi[1], bboxlo[2], xmin, ymin, zmin, xmax, ymax, zmax); - domain->lattice->bbox(1, bboxlo[0], bboxlo[1], bboxhi[2], xmin, ymin, zmin, xmax, ymax, zmax); - domain->lattice->bbox(1, bboxhi[0], bboxlo[1], bboxhi[2], xmin, ymin, zmin, xmax, ymax, zmax); - domain->lattice->bbox(1, bboxlo[0], bboxhi[1], bboxhi[2], xmin, ymin, zmin, xmax, ymax, zmax); - domain->lattice->bbox(1, bboxhi[0], bboxhi[1], bboxhi[2], xmin, ymin, zmin, xmax, ymax, zmax); + // for general triclinic, convert 8 corner points of bbox to general triclinic coords + // new set of 8 points is no longer an orthogonal bounding box + // instead invoke lattice->bbox() on each of 8 points + } else if (triclinic_general) { + double point[3]; + + point[0] = bboxlo[0]; point[1] = bboxlo[1]; point[2] = bboxlo[2]; + domain->restricted_to_general(point); + domain->lattice->bbox(1, point[0], point[1], point[2], xmin, ymin, zmin, xmax, ymax, zmax); + point[0] = bboxhi[0]; point[1] = bboxlo[1]; point[2] = bboxlo[2]; + domain->restricted_to_general(point); + domain->lattice->bbox(1, point[0], point[1], point[2], xmin, ymin, zmin, xmax, ymax, zmax); + + point[0] = bboxlo[0]; point[1] = bboxhi[1]; point[2] = bboxlo[2]; + domain->restricted_to_general(point); + domain->lattice->bbox(1, point[0], point[1], point[2], xmin, ymin, zmin, xmax, ymax, zmax); + point[0] = bboxhi[0]; point[1] = bboxhi[1]; point[2] = bboxlo[2]; + domain->restricted_to_general(point); + domain->lattice->bbox(1, point[0], point[1], point[2], xmin, ymin, zmin, xmax, ymax, zmax); + + point[0] = bboxlo[0]; point[1] = bboxlo[1]; point[2] = bboxhi[2]; + domain->restricted_to_general(point); + domain->lattice->bbox(1, point[0], point[1], point[2], xmin, ymin, zmin, xmax, ymax, zmax); + point[0] = bboxhi[0]; point[1] = bboxlo[1]; point[2] = bboxhi[2]; + domain->restricted_to_general(point); + domain->lattice->bbox(1, point[0], point[1], point[2], xmin, ymin, zmin, xmax, ymax, zmax); + + point[0] = bboxlo[0]; point[1] = bboxhi[1]; point[2] = bboxhi[2]; + domain->restricted_to_general(point); + domain->lattice->bbox(1, point[0], point[1], point[2], xmin, ymin, zmin, xmax, ymax, zmax); + point[0] = bboxhi[0]; point[1] = bboxhi[1]; point[2] = bboxhi[2]; + domain->restricted_to_general(point); + domain->lattice->bbox(1, point[0], point[1], point[2], xmin, ymin, zmin, xmax, ymax, zmax); + } + // ilo:ihi,jlo:jhi,klo:khi = loop bounds for lattice overlap of my subbox // overlap = any part of a unit cell (face,edge,pt) in common with my subbox // in lattice space, subbox is a tilted box @@ -1233,6 +1294,8 @@ void CreateAtoms::add_lattice() if (ymin < 0.0) jlo--; if (zmin < 0.0) klo--; + printf("LOOP LATTICE bounds: %d %d: %d %d: %d %d\n",ilo,ihi,jlo,jhi,klo,khi); + // count lattice sites on each proc nlatt_overflow = 0; @@ -1298,6 +1361,7 @@ void CreateAtoms::loop_lattice(int action) { int i, j, k, m; + int triclinic_general = domain->triclinic_general; const double *const *const basis = domain->lattice->basis; nlatt = 0; @@ -1317,6 +1381,10 @@ void CreateAtoms::loop_lattice(int action) domain->lattice->lattice2box(x[0], x[1], x[2]); + // convert from general to restricted triclinic coords + + if (triclinic_general) domain->general_to_restricted(x); + // if a region was specified, test if atom is in it if (style == REGION) diff --git a/src/create_atoms.h b/src/create_atoms.h index ae6f1b9d33..5850917112 100644 --- a/src/create_atoms.h +++ b/src/create_atoms.h @@ -69,6 +69,7 @@ class CreateAtoms : public Command { void add_mesh(const char *); int add_bisection(const double[3][3], tagint); int add_quasirandom(const double[3][3], tagint); + void add_lattice_triclinic_general(); void add_lattice(); void loop_lattice(int); void add_molecule(double *); diff --git a/src/create_box.cpp b/src/create_box.cpp index 74a8db8bb3..24dff6c01f 100644 --- a/src/create_box.cpp +++ b/src/create_box.cpp @@ -115,8 +115,10 @@ void CreateBox::command(int narg, char **arg) double chi = utils::numeric(FLERR, arg[iarg + 5], false, lmp); iarg += 6; - // use lattice2box() to generate avec, bvec, cvec and origin - + // use lattice2box() to generate origin and ABC vectors + // origin = abc lo + // ABC vector = hi in one dim - orign + double avec[3],bvec[3],cvec[3],origin[3]; double px,py,pz; @@ -144,12 +146,7 @@ void CreateBox::command(int narg, char **arg) 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 + // setup general triclinic box within Domain domain->set_general_triclinic(avec,bvec,cvec,origin); } diff --git a/src/domain.cpp b/src/domain.cpp index 9629645aa6..d86ced496b 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -667,7 +667,7 @@ void Domain::set_general_triclinic(double *avec_caller, double *bvec_caller, } /* ---------------------------------------------------------------------- - transform one atom's coords from general triclinic to restricted triclinic + transform atom coords from general triclinic to restricted triclinic ------------------------------------------------------------------------- */ void Domain::general_to_restricted(double *x) @@ -681,7 +681,7 @@ void Domain::general_to_restricted(double *x) } /* ---------------------------------------------------------------------- - transform one atom's coords from restricted triclinic to general triclinic + transform atom coords from restricted triclinic to general triclinic ------------------------------------------------------------------------- */ void Domain::restricted_to_general(double *x) diff --git a/src/lattice.cpp b/src/lattice.cpp index 539ec73531..2e7e740fb8 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -582,8 +582,9 @@ void Lattice::cross(double *x, double *y, double *z) } /* ---------------------------------------------------------------------- - convert x,y,z from lattice coords to box coords (flag = 0) or vice versa - use new point to expand bounding box (min to max) + convert x,y,z from lattice coords to box coords (flag = 0) + or from box coords to lattice coords (flag = 1) + either way, use new point to expand bounding box (min to max) ------------------------------------------------------------------------- */ void Lattice::bbox(int flag, double x, double y, double z,