diff --git a/src/domain.cpp b/src/domain.cpp
index 0590b956d..2c3f61401 100644
--- a/src/domain.cpp
+++ b/src/domain.cpp
@@ -1,1930 +1,1933 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author (triclinic) : Pieter in 't Veld (SNL)
------------------------------------------------------------------------- */
#include <mpi.h>
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <math.h>
#include "domain.h"
#include "style_region.h"
#include "atom.h"
#include "atom_vec.h"
#include "molecule.h"
#include "force.h"
#include "kspace.h"
#include "update.h"
#include "modify.h"
#include "fix.h"
#include "fix_deform.h"
#include "region.h"
#include "lattice.h"
#include "comm.h"
#include "output.h"
#include "thermo.h"
#include "universe.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace MathConst;
enum{NO_REMAP,X_REMAP,V_REMAP}; // same as fix_deform.cpp
enum{IGNORE,WARN,ERROR}; // same as thermo.cpp
enum{LAYOUT_UNIFORM,LAYOUT_NONUNIFORM,LAYOUT_TILED}; // several files
#define BIG 1.0e20
#define SMALL 1.0e-4
#define DELTAREGION 4
#define BONDSTRETCH 1.1
/* ----------------------------------------------------------------------
default is periodic
------------------------------------------------------------------------- */
Domain::Domain(LAMMPS *lmp) : Pointers(lmp)
{
box_exist = 0;
dimension = 3;
nonperiodic = 0;
xperiodic = yperiodic = zperiodic = 1;
periodicity[0] = xperiodic;
periodicity[1] = yperiodic;
periodicity[2] = zperiodic;
boundary[0][0] = boundary[0][1] = 0;
boundary[1][0] = boundary[1][1] = 0;
boundary[2][0] = boundary[2][1] = 0;
minxlo = minxhi = 0.0;
minylo = minyhi = 0.0;
minzlo = minzhi = 0.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;
h[3] = h[4] = h[5] = 0.0;
h_inv[3] = h_inv[4] = h_inv[5] = 0.0;
h_rate[0] = h_rate[1] = h_rate[2] =
h_rate[3] = h_rate[4] = h_rate[5] = 0.0;
h_ratelo[0] = h_ratelo[1] = h_ratelo[2] = 0.0;
prd_lamda[0] = prd_lamda[1] = prd_lamda[2] = 1.0;
prd_half_lamda[0] = prd_half_lamda[1] = prd_half_lamda[2] = 0.5;
boxlo_lamda[0] = boxlo_lamda[1] = boxlo_lamda[2] = 0.0;
boxhi_lamda[0] = boxhi_lamda[1] = boxhi_lamda[2] = 1.0;
lattice = NULL;
char **args = new char*[2];
args[0] = (char *) "none";
args[1] = (char *) "1.0";
set_lattice(2,args);
delete [] args;
nregion = maxregion = 0;
regions = NULL;
copymode = 0;
+
+ region_map = new RegionCreatorMap();
+
+#define REGION_CLASS
+#define RegionStyle(key,Class) \
+ (*region_map)[#key] = ®ion_creator<Class>;
+#include "style_region.h"
+#undef RegionStyle
+#undef REGION_CLASS
}
/* ---------------------------------------------------------------------- */
Domain::~Domain()
{
if (copymode) return;
delete lattice;
for (int i = 0; i < nregion; i++) delete regions[i];
memory->sfree(regions);
+
+ delete region_map;
}
/* ---------------------------------------------------------------------- */
void Domain::init()
{
// set box_change flags if box size/shape/sub-domains ever change
// due to shrink-wrapping or fixes that change box size/shape/sub-domains
box_change_size = box_change_shape = box_change_domain = 0;
if (nonperiodic == 2) box_change_size = 1;
for (int i = 0; i < modify->nfix; i++) {
if (modify->fix[i]->box_change_size) box_change_size = 1;
if (modify->fix[i]->box_change_shape) box_change_shape = 1;
if (modify->fix[i]->box_change_domain) box_change_domain = 1;
}
box_change = 0;
if (box_change_size || box_change_shape || box_change_domain) box_change = 1;
// check for fix deform
deform_flag = deform_vremap = deform_groupbit = 0;
for (int i = 0; i < modify->nfix; i++)
if (strcmp(modify->fix[i]->style,"deform") == 0) {
deform_flag = 1;
if (((FixDeform *) modify->fix[i])->remapflag == V_REMAP) {
deform_vremap = 1;
deform_groupbit = modify->fix[i]->groupbit;
}
}
// region inits
for (int i = 0; i < nregion; i++) regions[i]->init();
}
/* ----------------------------------------------------------------------
set initial global box
assumes boxlo/hi and triclinic tilts are already set
expandflag = 1 if need to expand box in shrink-wrapped dims
not invoked by read_restart since box is already expanded
if don't prevent further expansion, restarted triclinic box
with unchanged tilt factors can become a box with atoms outside the box
------------------------------------------------------------------------- */
void Domain::set_initial_box(int expandflag)
{
// error checks for orthogonal and triclinic domains
if (boxlo[0] >= boxhi[0] || boxlo[1] >= boxhi[1] || boxlo[2] >= boxhi[2])
error->one(FLERR,"Box bounds are invalid or missing");
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 ((fabs(xy/(boxhi[0]-boxlo[0])) > 0.5 && xperiodic) ||
(fabs(xz/(boxhi[0]-boxlo[0])) > 0.5 && xperiodic) ||
(fabs(yz/(boxhi[1]-boxlo[1])) > 0.5 && yperiodic)) {
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
// this works for any unit system
small[0] = SMALL * (boxhi[0] - boxlo[0]);
small[1] = SMALL * (boxhi[1] - boxlo[1]);
small[2] = SMALL * (boxhi[2] - boxlo[2]);
// if expandflag, adjust box lo/hi for shrink-wrapped dims
if (!expandflag) return;
if (boundary[0][0] == 2) boxlo[0] -= small[0];
else if (boundary[0][0] == 3) minxlo = boxlo[0];
if (boundary[0][1] == 2) boxhi[0] += small[0];
else if (boundary[0][1] == 3) minxhi = boxhi[0];
if (boundary[1][0] == 2) boxlo[1] -= small[1];
else if (boundary[1][0] == 3) minylo = boxlo[1];
if (boundary[1][1] == 2) boxhi[1] += small[1];
else if (boundary[1][1] == 3) minyhi = boxhi[1];
if (boundary[2][0] == 2) boxlo[2] -= small[2];
else if (boundary[2][0] == 3) minzlo = boxlo[2];
if (boundary[2][1] == 2) boxhi[2] += small[2];
else if (boundary[2][1] == 3) minzhi = boxhi[2];
}
/* ----------------------------------------------------------------------
set global box params
assumes boxlo/hi and triclinic tilts are already set
------------------------------------------------------------------------- */
void Domain::set_global_box()
{
prd[0] = xprd = boxhi[0] - boxlo[0];
prd[1] = yprd = boxhi[1] - boxlo[1];
prd[2] = zprd = boxhi[2] - boxlo[2];
h[0] = xprd;
h[1] = yprd;
h[2] = zprd;
h_inv[0] = 1.0/h[0];
h_inv[1] = 1.0/h[1];
h_inv[2] = 1.0/h[2];
prd_half[0] = xprd_half = 0.5*xprd;
prd_half[1] = yprd_half = 0.5*yprd;
prd_half[2] = zprd_half = 0.5*zprd;
if (triclinic) {
h[3] = yz;
h[4] = xz;
h[5] = xy;
h_inv[3] = -h[3] / (h[1]*h[2]);
h_inv[4] = (h[3]*h[5] - h[1]*h[4]) / (h[0]*h[1]*h[2]);
h_inv[5] = -h[5] / (h[0]*h[1]);
boxlo_bound[0] = MIN(boxlo[0],boxlo[0]+xy);
boxlo_bound[0] = MIN(boxlo_bound[0],boxlo_bound[0]+xz);
boxlo_bound[1] = MIN(boxlo[1],boxlo[1]+yz);
boxlo_bound[2] = boxlo[2];
boxhi_bound[0] = MAX(boxhi[0],boxhi[0]+xy);
boxhi_bound[0] = MAX(boxhi_bound[0],boxhi_bound[0]+xz);
boxhi_bound[1] = MAX(boxhi[1],boxhi[1]+yz);
boxhi_bound[2] = boxhi[2];
}
}
/* ----------------------------------------------------------------------
set lamda box params
assumes global box is defined and proc assignment has been made
uses comm->xyz_split or comm->mysplit
to define subbox boundaries in consistent manner
------------------------------------------------------------------------- */
void Domain::set_lamda_box()
{
if (comm->layout != LAYOUT_TILED) {
int *myloc = comm->myloc;
double *xsplit = comm->xsplit;
double *ysplit = comm->ysplit;
double *zsplit = comm->zsplit;
sublo_lamda[0] = xsplit[myloc[0]];
subhi_lamda[0] = xsplit[myloc[0]+1];
sublo_lamda[1] = ysplit[myloc[1]];
subhi_lamda[1] = ysplit[myloc[1]+1];
sublo_lamda[2] = zsplit[myloc[2]];
subhi_lamda[2] = zsplit[myloc[2]+1];
} else {
double (*mysplit)[2] = comm->mysplit;
sublo_lamda[0] = mysplit[0][0];
subhi_lamda[0] = mysplit[0][1];
sublo_lamda[1] = mysplit[1][0];
subhi_lamda[1] = mysplit[1][1];
sublo_lamda[2] = mysplit[2][0];
subhi_lamda[2] = mysplit[2][1];
}
}
/* ----------------------------------------------------------------------
set local subbox params for orthogonal boxes
assumes global box is defined and proc assignment has been made
uses comm->xyz_split or comm->mysplit
to define subbox boundaries in consistent manner
insure subhi[max] = boxhi
------------------------------------------------------------------------- */
void Domain::set_local_box()
{
if (triclinic) return;
if (comm->layout != LAYOUT_TILED) {
int *myloc = comm->myloc;
int *procgrid = comm->procgrid;
double *xsplit = comm->xsplit;
double *ysplit = comm->ysplit;
double *zsplit = comm->zsplit;
sublo[0] = boxlo[0] + xprd*xsplit[myloc[0]];
if (myloc[0] < procgrid[0]-1) subhi[0] = boxlo[0] + xprd*xsplit[myloc[0]+1];
else subhi[0] = boxhi[0];
sublo[1] = boxlo[1] + yprd*ysplit[myloc[1]];
if (myloc[1] < procgrid[1]-1) subhi[1] = boxlo[1] + yprd*ysplit[myloc[1]+1];
else subhi[1] = boxhi[1];
sublo[2] = boxlo[2] + zprd*zsplit[myloc[2]];
if (myloc[2] < procgrid[2]-1) subhi[2] = boxlo[2] + zprd*zsplit[myloc[2]+1];
else subhi[2] = boxhi[2];
} else {
double (*mysplit)[2] = comm->mysplit;
sublo[0] = boxlo[0] + xprd*mysplit[0][0];
if (mysplit[0][1] < 1.0) subhi[0] = boxlo[0] + xprd*mysplit[0][1];
else subhi[0] = boxhi[0];
sublo[1] = boxlo[1] + yprd*mysplit[1][0];
if (mysplit[1][1] < 1.0) subhi[1] = boxlo[1] + yprd*mysplit[1][1];
else subhi[1] = boxhi[1];
sublo[2] = boxlo[2] + zprd*mysplit[2][0];
if (mysplit[2][1] < 1.0) subhi[2] = boxlo[2] + zprd*mysplit[2][1];
else subhi[2] = boxhi[2];
}
}
/* ----------------------------------------------------------------------
reset global & local boxes due to global box boundary changes
if shrink-wrapped, determine atom extent and reset boxlo/hi
for triclinic, atoms must be in lamda coords (0-1) before reset_box is called
------------------------------------------------------------------------- */
void Domain::reset_box()
{
// perform shrink-wrapping
// compute extent of atoms on this proc
// for triclinic, this is done in lamda space
if (nonperiodic == 2) {
double extent[3][2],all[3][2];
extent[2][0] = extent[1][0] = extent[0][0] = BIG;
extent[2][1] = extent[1][1] = extent[0][1] = -BIG;
double **x = atom->x;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
extent[0][0] = MIN(extent[0][0],x[i][0]);
extent[0][1] = MAX(extent[0][1],x[i][0]);
extent[1][0] = MIN(extent[1][0],x[i][1]);
extent[1][1] = MAX(extent[1][1],x[i][1]);
extent[2][0] = MIN(extent[2][0],x[i][2]);
extent[2][1] = MAX(extent[2][1],x[i][2]);
}
// compute extent across all procs
// flip sign of MIN to do it in one Allreduce MAX
extent[0][0] = -extent[0][0];
extent[1][0] = -extent[1][0];
extent[2][0] = -extent[2][0];
MPI_Allreduce(extent,all,6,MPI_DOUBLE,MPI_MAX,world);
// for triclinic, convert back to box coords before changing box
if (triclinic) lamda2x(atom->nlocal);
// in shrink-wrapped dims, set box by atom extent
// if minimum set, enforce min box size settings
// for triclinic, convert lamda extent to box coords, then set box lo/hi
// decided NOT to do the next comment - don't want to sneakily change tilt
// for triclinic, adjust tilt factors if 2nd dim is shrink-wrapped,
// so that displacement in 1st dim stays the same
if (triclinic == 0) {
if (xperiodic == 0) {
if (boundary[0][0] == 2) boxlo[0] = -all[0][0] - small[0];
else if (boundary[0][0] == 3)
boxlo[0] = MIN(-all[0][0]-small[0],minxlo);
if (boundary[0][1] == 2) boxhi[0] = all[0][1] + small[0];
else if (boundary[0][1] == 3) boxhi[0] = MAX(all[0][1]+small[0],minxhi);
if (boxlo[0] > boxhi[0]) error->all(FLERR,"Illegal simulation box");
}
if (yperiodic == 0) {
if (boundary[1][0] == 2) boxlo[1] = -all[1][0] - small[1];
else if (boundary[1][0] == 3)
boxlo[1] = MIN(-all[1][0]-small[1],minylo);
if (boundary[1][1] == 2) boxhi[1] = all[1][1] + small[1];
else if (boundary[1][1] == 3) boxhi[1] = MAX(all[1][1]+small[1],minyhi);
if (boxlo[1] > boxhi[1]) error->all(FLERR,"Illegal simulation box");
}
if (zperiodic == 0) {
if (boundary[2][0] == 2) boxlo[2] = -all[2][0] - small[2];
else if (boundary[2][0] == 3)
boxlo[2] = MIN(-all[2][0]-small[2],minzlo);
if (boundary[2][1] == 2) boxhi[2] = all[2][1] + small[2];
else if (boundary[2][1] == 3) boxhi[2] = MAX(all[2][1]+small[2],minzhi);
if (boxlo[2] > boxhi[2]) error->all(FLERR,"Illegal simulation box");
}
} else {
double lo[3],hi[3];
if (xperiodic == 0) {
lo[0] = -all[0][0]; lo[1] = 0.0; lo[2] = 0.0;
lamda2x(lo,lo);
hi[0] = all[0][1]; hi[1] = 0.0; hi[2] = 0.0;
lamda2x(hi,hi);
if (boundary[0][0] == 2) boxlo[0] = lo[0] - small[0];
else if (boundary[0][0] == 3) boxlo[0] = MIN(lo[0]-small[0],minxlo);
if (boundary[0][1] == 2) boxhi[0] = hi[0] + small[0];
else if (boundary[0][1] == 3) boxhi[0] = MAX(hi[0]+small[0],minxhi);
if (boxlo[0] > boxhi[0]) error->all(FLERR,"Illegal simulation box");
}
if (yperiodic == 0) {
lo[0] = 0.0; lo[1] = -all[1][0]; lo[2] = 0.0;
lamda2x(lo,lo);
hi[0] = 0.0; hi[1] = all[1][1]; hi[2] = 0.0;
lamda2x(hi,hi);
if (boundary[1][0] == 2) boxlo[1] = lo[1] - small[1];
else if (boundary[1][0] == 3) boxlo[1] = MIN(lo[1]-small[1],minylo);
if (boundary[1][1] == 2) boxhi[1] = hi[1] + small[1];
else if (boundary[1][1] == 3) boxhi[1] = MAX(hi[1]+small[1],minyhi);
if (boxlo[1] > boxhi[1]) error->all(FLERR,"Illegal simulation box");
//xy *= (boxhi[1]-boxlo[1]) / yprd;
}
if (zperiodic == 0) {
lo[0] = 0.0; lo[1] = 0.0; lo[2] = -all[2][0];
lamda2x(lo,lo);
hi[0] = 0.0; hi[1] = 0.0; hi[2] = all[2][1];
lamda2x(hi,hi);
if (boundary[2][0] == 2) boxlo[2] = lo[2] - small[2];
else if (boundary[2][0] == 3) boxlo[2] = MIN(lo[2]-small[2],minzlo);
if (boundary[2][1] == 2) boxhi[2] = hi[2] + small[2];
else if (boundary[2][1] == 3) boxhi[2] = MAX(hi[2]+small[2],minzhi);
if (boxlo[2] > boxhi[2]) error->all(FLERR,"Illegal simulation box");
//xz *= (boxhi[2]-boxlo[2]) / xprd;
//yz *= (boxhi[2]-boxlo[2]) / yprd;
}
}
}
// reset box whether shrink-wrapping or not
set_global_box();
set_local_box();
// if shrink-wrapped & kspace is defined (i.e. using MSM), call setup()
// also call init() (to test for compatibility) ?
if (nonperiodic == 2 && force->kspace) {
//force->kspace->init();
force->kspace->setup();
}
// if shrink-wrapped & triclinic, re-convert to lamda coords for new box
// re-invoke pbc() b/c x2lamda result can be outside [0,1] due to roundoff
if (nonperiodic == 2 && triclinic) {
x2lamda(atom->nlocal);
pbc();
}
}
/* ----------------------------------------------------------------------
enforce PBC and modify box image flags for each atom
called every reneighboring and by other commands that change atoms
resulting coord must satisfy lo <= coord < hi
MAX is important since coord - prd < lo can happen when coord = hi
if fix deform, remap velocity of fix group atoms by box edge velocities
for triclinic, atoms must be in lamda coords (0-1) before pbc is called
image = 10 or 20 bits for each dimension depending on sizeof(imageint)
increment/decrement in wrap-around fashion
------------------------------------------------------------------------- */
void Domain::pbc()
{
int i;
imageint idim,otherdims;
double *lo,*hi,*period;
int nlocal = atom->nlocal;
double **x = atom->x;
double **v = atom->v;
int *mask = atom->mask;
imageint *image = atom->image;
// verify owned atoms have valid numerical coords
// may not if computed pairwise force between 2 atoms at same location
double *coord;
int n3 = 3*nlocal;
coord = &x[0][0]; // note: x is always initialzed to at least one element.
int flag = 0;
for (i = 0; i < n3; i++)
if (!ISFINITE(*coord++)) flag = 1;
if (flag) error->one(FLERR,"Non-numeric atom coords - simulation unstable");
// setup for PBC checks
if (triclinic == 0) {
lo = boxlo;
hi = boxhi;
period = prd;
} else {
lo = boxlo_lamda;
hi = boxhi_lamda;
period = prd_lamda;
}
// apply PBC to each owned atom
for (i = 0; i < nlocal; i++) {
if (xperiodic) {
if (x[i][0] < lo[0]) {
x[i][0] += period[0];
if (deform_vremap && mask[i] & deform_groupbit) v[i][0] += h_rate[0];
idim = image[i] & IMGMASK;
otherdims = image[i] ^ idim;
idim--;
idim &= IMGMASK;
image[i] = otherdims | idim;
}
if (x[i][0] >= hi[0]) {
x[i][0] -= period[0];
x[i][0] = MAX(x[i][0],lo[0]);
if (deform_vremap && mask[i] & deform_groupbit) v[i][0] -= h_rate[0];
idim = image[i] & IMGMASK;
otherdims = image[i] ^ idim;
idim++;
idim &= IMGMASK;
image[i] = otherdims | idim;
}
}
if (yperiodic) {
if (x[i][1] < lo[1]) {
x[i][1] += period[1];
if (deform_vremap && mask[i] & deform_groupbit) {
v[i][0] += h_rate[5];
v[i][1] += h_rate[1];
}
idim = (image[i] >> IMGBITS) & IMGMASK;
otherdims = image[i] ^ (idim << IMGBITS);
idim--;
idim &= IMGMASK;
image[i] = otherdims | (idim << IMGBITS);
}
if (x[i][1] >= hi[1]) {
x[i][1] -= period[1];
x[i][1] = MAX(x[i][1],lo[1]);
if (deform_vremap && mask[i] & deform_groupbit) {
v[i][0] -= h_rate[5];
v[i][1] -= h_rate[1];
}
idim = (image[i] >> IMGBITS) & IMGMASK;
otherdims = image[i] ^ (idim << IMGBITS);
idim++;
idim &= IMGMASK;
image[i] = otherdims | (idim << IMGBITS);
}
}
if (zperiodic) {
if (x[i][2] < lo[2]) {
x[i][2] += period[2];
if (deform_vremap && mask[i] & deform_groupbit) {
v[i][0] += h_rate[4];
v[i][1] += h_rate[3];
v[i][2] += h_rate[2];
}
idim = image[i] >> IMG2BITS;
otherdims = image[i] ^ (idim << IMG2BITS);
idim--;
idim &= IMGMASK;
image[i] = otherdims | (idim << IMG2BITS);
}
if (x[i][2] >= hi[2]) {
x[i][2] -= period[2];
x[i][2] = MAX(x[i][2],lo[2]);
if (deform_vremap && mask[i] & deform_groupbit) {
v[i][0] -= h_rate[4];
v[i][1] -= h_rate[3];
v[i][2] -= h_rate[2];
}
idim = image[i] >> IMG2BITS;
otherdims = image[i] ^ (idim << IMG2BITS);
idim++;
idim &= IMGMASK;
image[i] = otherdims | (idim << IMG2BITS);
}
}
}
}
/* ----------------------------------------------------------------------
check that point is inside box boundaries, in [lo,hi) sense
return 1 if true, 0 if false
------------------------------------------------------------------------- */
int Domain::inside(double* x)
{
double *lo,*hi;
double lamda[3];
if (triclinic == 0) {
lo = boxlo;
hi = boxhi;
if (x[0] < lo[0] || x[0] >= hi[0] ||
x[1] < lo[1] || x[1] >= hi[1] ||
x[2] < lo[2] || x[2] >= hi[2]) return 0;
else return 1;
} else {
lo = boxlo_lamda;
hi = boxhi_lamda;
x2lamda(x,lamda);
if (lamda[0] < lo[0] || lamda[0] >= hi[0] ||
lamda[1] < lo[1] || lamda[1] >= hi[1] ||
lamda[2] < lo[2] || lamda[2] >= hi[2]) return 0;
else return 1;
}
}
/* ----------------------------------------------------------------------
check that point is inside nonperiodic boundaries, in [lo,hi) sense
return 1 if true, 0 if false
------------------------------------------------------------------------- */
int Domain::inside_nonperiodic(double* x)
{
double *lo,*hi;
double lamda[3];
if (xperiodic && yperiodic && zperiodic) return 1;
if (triclinic == 0) {
lo = boxlo;
hi = boxhi;
if (!xperiodic && (x[0] < lo[0] || x[0] >= hi[0])) return 0;
if (!yperiodic && (x[1] < lo[1] || x[1] >= hi[1])) return 0;
if (!zperiodic && (x[2] < lo[2] || x[2] >= hi[2])) return 0;
return 1;
} else {
lo = boxlo_lamda;
hi = boxhi_lamda;
x2lamda(x,lamda);
if (!xperiodic && (lamda[0] < lo[0] || lamda[0] >= hi[0])) return 0;
if (!yperiodic && (lamda[1] < lo[1] || lamda[1] >= hi[1])) return 0;
if (!zperiodic && (lamda[2] < lo[2] || lamda[2] >= hi[2])) return 0;
return 1;
}
}
/* ----------------------------------------------------------------------
warn if image flags of any bonded atoms are inconsistent
could be a problem when using replicate or fix rigid
------------------------------------------------------------------------- */
void Domain::image_check()
{
int i,j,k,n,imol,iatom;
tagint tagprev;
// only need to check if system is molecular and some dimension is periodic
// if running verlet/split, don't check on KSpace partition since
// it has no ghost atoms and thus bond partners won't exist
if (!atom->molecular) return;
if (!xperiodic && !yperiodic && (dimension == 2 || !zperiodic)) return;
if (strncmp(update->integrate_style,"verlet/split",12) == 0 &&
universe->iworld != 0) return;
// communicate unwrapped position of owned atoms to ghost atoms
double **unwrap;
memory->create(unwrap,atom->nmax,3,"domain:unwrap");
double **x = atom->x;
imageint *image = atom->image;
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++)
unmap(x[i],image[i],unwrap[i]);
comm->forward_comm_array(3,unwrap);
// compute unwrapped extent of each bond
// flag if any bond component is longer than 1/2 of periodic box length
// flag if any bond component is longer than non-periodic box length
// which means image flags in that dimension were different
int molecular = atom->molecular;
int *num_bond = atom->num_bond;
tagint **bond_atom = atom->bond_atom;
int **bond_type = atom->bond_type;
tagint *tag = atom->tag;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
Molecule **onemols = atom->avec->onemols;
double delx,dely,delz;
int lostbond = output->thermo->lostbond;
int nmissing = 0;
int flag = 0;
for (i = 0; i < nlocal; i++) {
if (molecular == 1) n = num_bond[i];
else {
if (molindex[i] < 0) continue;
imol = molindex[i];
iatom = molatom[i];
n = onemols[imol]->num_bond[iatom];
}
for (j = 0; j < n; j++) {
if (molecular == 1) {
if (bond_type[i][j] <= 0) continue;
k = atom->map(bond_atom[i][j]);
} else {
if (onemols[imol]->bond_type[iatom][j] < 0) continue;
tagprev = tag[i] - iatom - 1;
k = atom->map(onemols[imol]->bond_atom[iatom][j]+tagprev);
}
if (k == -1) {
nmissing++;
if (lostbond == ERROR)
error->one(FLERR,"Bond atom missing in image check");
continue;
}
delx = unwrap[i][0] - unwrap[k][0];
dely = unwrap[i][1] - unwrap[k][1];
delz = unwrap[i][2] - unwrap[k][2];
if (xperiodic && delx > xprd_half) flag = 1;
if (xperiodic && dely > yprd_half) flag = 1;
if (dimension == 3 && zperiodic && delz > zprd_half) flag = 1;
if (!xperiodic && delx > xprd) flag = 1;
if (!yperiodic && dely > yprd) flag = 1;
if (dimension == 3 && !zperiodic && delz > zprd) flag = 1;
}
}
int flagall;
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_MAX,world);
if (flagall && comm->me == 0)
error->warning(FLERR,"Inconsistent image flags");
if (lostbond == WARN) {
int all;
MPI_Allreduce(&nmissing,&all,1,MPI_INT,MPI_SUM,world);
if (all && comm->me == 0)
error->warning(FLERR,"Bond atom missing in image check");
}
memory->destroy(unwrap);
}
/* ----------------------------------------------------------------------
warn if end atoms in any bonded interaction
are further apart than half a periodic box length
could cause problems when bonded neighbor list is built since
closest_image() could return wrong image
------------------------------------------------------------------------- */
void Domain::box_too_small_check()
{
int i,j,k,n,imol,iatom;
tagint tagprev;
// only need to check if system is molecular and some dimension is periodic
// if running verlet/split, don't check on KSpace partition since
// it has no ghost atoms and thus bond partners won't exist
if (!atom->molecular) return;
if (!xperiodic && !yperiodic && (dimension == 2 || !zperiodic)) return;
if (strncmp(update->integrate_style,"verlet/split",12) == 0 &&
universe->iworld != 0) return;
// maxbondall = longest current bond length
// if periodic box dim is tiny (less than 2 * bond-length),
// minimum_image() itself may compute bad bond lengths
// in this case, image_check() should warn,
// assuming 2 atoms have consistent image flags
int molecular = atom->molecular;
double **x = atom->x;
int *num_bond = atom->num_bond;
tagint **bond_atom = atom->bond_atom;
int **bond_type = atom->bond_type;
tagint *tag = atom->tag;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
Molecule **onemols = atom->avec->onemols;
int nlocal = atom->nlocal;
double delx,dely,delz,rsq;
double maxbondme = 0.0;
int lostbond = output->thermo->lostbond;
int nmissing = 0;
for (i = 0; i < nlocal; i++) {
if (molecular == 1) n = num_bond[i];
else {
if (molindex[i] < 0) continue;
imol = molindex[i];
iatom = molatom[i];
n = onemols[imol]->num_bond[iatom];
}
for (j = 0; j < n; j++) {
if (molecular == 1) {
if (bond_type[i][j] <= 0) continue;
k = atom->map(bond_atom[i][j]);
} else {
if (onemols[imol]->bond_type[iatom][j] < 0) continue;
tagprev = tag[i] - iatom - 1;
k = atom->map(onemols[imol]->bond_atom[iatom][j]+tagprev);
}
if (k == -1) {
nmissing++;
if (lostbond == ERROR)
error->one(FLERR,"Bond atom missing in box size check");
continue;
}
delx = x[i][0] - x[k][0];
dely = x[i][1] - x[k][1];
delz = x[i][2] - x[k][2];
minimum_image(delx,dely,delz);
rsq = delx*delx + dely*dely + delz*delz;
maxbondme = MAX(maxbondme,rsq);
}
}
if (lostbond == WARN) {
int all;
MPI_Allreduce(&nmissing,&all,1,MPI_INT,MPI_SUM,world);
if (all && comm->me == 0)
error->warning(FLERR,"Bond atom missing in box size check");
}
double maxbondall;
MPI_Allreduce(&maxbondme,&maxbondall,1,MPI_DOUBLE,MPI_MAX,world);
maxbondall = sqrt(maxbondall);
// maxdelta = furthest apart 2 atoms in a bonded interaction can be
// include BONDSTRETCH factor to account for dynamics
double maxdelta = maxbondall * BONDSTRETCH;
if (atom->nangles) maxdelta = 2.0 * maxbondall * BONDSTRETCH;
if (atom->ndihedrals) maxdelta = 3.0 * maxbondall * BONDSTRETCH;
// warn if maxdelta > than half any periodic box length
// since atoms in the interaction could rotate into that dimension
int flag = 0;
if (xperiodic && maxdelta > xprd_half) flag = 1;
if (yperiodic && maxdelta > yprd_half) flag = 1;
if (dimension == 3 && zperiodic && maxdelta > zprd_half) flag = 1;
int flagall;
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_MAX,world);
if (flagall && comm->me == 0)
error->warning(FLERR,
"Bond/angle/dihedral extent > half of periodic box length");
}
/* ----------------------------------------------------------------------
check warn if any proc's subbox is smaller than thresh
since may lead to lost atoms in comm->exchange()
current callers set thresh = neighbor skin
------------------------------------------------------------------------- */
void Domain::subbox_too_small_check(double thresh)
{
int flag = 0;
if (!triclinic) {
if (subhi[0]-sublo[0] < thresh || subhi[1]-sublo[1] < thresh) flag = 1;
if (dimension == 3 && subhi[2]-sublo[2] < thresh) flag = 1;
} else {
double delta = subhi_lamda[0] - sublo_lamda[0];
if (delta*prd[0] < thresh) flag = 1;
delta = subhi_lamda[1] - sublo_lamda[1];
if (delta*prd[1] < thresh) flag = 1;
if (dimension == 3) {
delta = subhi_lamda[2] - sublo_lamda[2];
if (delta*prd[2] < thresh) flag = 1;
}
}
int flagall;
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
if (flagall && comm->me == 0)
error->warning(FLERR,"Proc sub-domain size < neighbor skin, "
"could lead to lost atoms");
}
/* ----------------------------------------------------------------------
minimum image convention
use 1/2 of box size as test
for triclinic, also add/subtract tilt factors in other dims as needed
------------------------------------------------------------------------- */
void Domain::minimum_image(double &dx, double &dy, double &dz)
{
if (triclinic == 0) {
if (xperiodic) {
if (fabs(dx) > xprd_half) {
if (dx < 0.0) dx += xprd;
else dx -= xprd;
}
}
if (yperiodic) {
if (fabs(dy) > yprd_half) {
if (dy < 0.0) dy += yprd;
else dy -= yprd;
}
}
if (zperiodic) {
if (fabs(dz) > zprd_half) {
if (dz < 0.0) dz += zprd;
else dz -= zprd;
}
}
} else {
if (zperiodic) {
if (fabs(dz) > zprd_half) {
if (dz < 0.0) {
dz += zprd;
dy += yz;
dx += xz;
} else {
dz -= zprd;
dy -= yz;
dx -= xz;
}
}
}
if (yperiodic) {
if (fabs(dy) > yprd_half) {
if (dy < 0.0) {
dy += yprd;
dx += xy;
} else {
dy -= yprd;
dx -= xy;
}
}
}
if (xperiodic) {
if (fabs(dx) > xprd_half) {
if (dx < 0.0) dx += xprd;
else dx -= xprd;
}
}
}
}
/* ----------------------------------------------------------------------
minimum image convention
use 1/2 of box size as test
for triclinic, also add/subtract tilt factors in other dims as needed
------------------------------------------------------------------------- */
void Domain::minimum_image(double *delta)
{
if (triclinic == 0) {
if (xperiodic) {
if (fabs(delta[0]) > xprd_half) {
if (delta[0] < 0.0) delta[0] += xprd;
else delta[0] -= xprd;
}
}
if (yperiodic) {
if (fabs(delta[1]) > yprd_half) {
if (delta[1] < 0.0) delta[1] += yprd;
else delta[1] -= yprd;
}
}
if (zperiodic) {
if (fabs(delta[2]) > zprd_half) {
if (delta[2] < 0.0) delta[2] += zprd;
else delta[2] -= zprd;
}
}
} else {
if (zperiodic) {
if (fabs(delta[2]) > zprd_half) {
if (delta[2] < 0.0) {
delta[2] += zprd;
delta[1] += yz;
delta[0] += xz;
} else {
delta[2] -= zprd;
delta[1] -= yz;
delta[0] -= xz;
}
}
}
if (yperiodic) {
if (fabs(delta[1]) > yprd_half) {
if (delta[1] < 0.0) {
delta[1] += yprd;
delta[0] += xy;
} else {
delta[1] -= yprd;
delta[0] -= xy;
}
}
}
if (xperiodic) {
if (fabs(delta[0]) > xprd_half) {
if (delta[0] < 0.0) delta[0] += xprd;
else delta[0] -= xprd;
}
}
}
}
/* ----------------------------------------------------------------------
return local index of atom J or any of its images that is closest to atom I
if J is not a valid index like -1, just return it
------------------------------------------------------------------------- */
int Domain::closest_image(int i, int j)
{
if (j < 0) return j;
int *sametag = atom->sametag;
double **x = atom->x;
double *xi = x[i];
int closest = j;
double delx = xi[0] - x[j][0];
double dely = xi[1] - x[j][1];
double delz = xi[2] - x[j][2];
double rsqmin = delx*delx + dely*dely + delz*delz;
double rsq;
while (sametag[j] >= 0) {
j = sametag[j];
delx = xi[0] - x[j][0];
dely = xi[1] - x[j][1];
delz = xi[2] - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq < rsqmin) {
rsqmin = rsq;
closest = j;
}
}
return closest;
}
/* ----------------------------------------------------------------------
find and return Xj image = periodic image of Xj that is closest to Xi
for triclinic, add/subtract tilt factors in other dims as needed
------------------------------------------------------------------------- */
void Domain::closest_image(const double * const xi, const double * const xj,
double * const xjimage)
{
double dx = xj[0] - xi[0];
double dy = xj[1] - xi[1];
double dz = xj[2] - xi[2];
if (triclinic == 0) {
if (xperiodic) {
if (dx < 0.0) {
while (dx < 0.0) dx += xprd;
if (dx > xprd_half) dx -= xprd;
} else {
while (dx > 0.0) dx -= xprd;
if (dx < -xprd_half) dx += xprd;
}
}
if (yperiodic) {
if (dy < 0.0) {
while (dy < 0.0) dy += yprd;
if (dy > yprd_half) dy -= yprd;
} else {
while (dy > 0.0) dy -= yprd;
if (dy < -yprd_half) dy += yprd;
}
}
if (zperiodic) {
if (dz < 0.0) {
while (dz < 0.0) dz += zprd;
if (dz > zprd_half) dz -= zprd;
} else {
while (dz > 0.0) dz -= zprd;
if (dz < -zprd_half) dz += zprd;
}
}
} else {
if (zperiodic) {
if (dz < 0.0) {
while (dz < 0.0) {
dz += zprd;
dy += yz;
dx += xz;
}
if (dz > zprd_half) {
dz -= zprd;
dy -= yz;
dx -= xz;
}
} else {
while (dz > 0.0) {
dz -= zprd;
dy -= yz;
dx -= xz;
}
if (dz < -zprd_half) {
dz += zprd;
dy += yz;
dx += xz;
}
}
}
if (yperiodic) {
if (dy < 0.0) {
while (dy < 0.0) {
dy += yprd;
dx += xy;
}
if (dy > yprd_half) {
dy -= yprd;
dx -= xy;
}
} else {
while (dy > 0.0) {
dy -= yprd;
dx -= xy;
}
if (dy < -yprd_half) {
dy += yprd;
dx += xy;
}
}
}
if (xperiodic) {
if (dx < 0.0) {
while (dx < 0.0) dx += xprd;
if (dx > xprd_half) dx -= xprd;
} else {
while (dx > 0.0) dx -= xprd;
if (dx < -xprd_half) dx += xprd;
}
}
}
xjimage[0] = xi[0] + dx;
xjimage[1] = xi[1] + dy;
xjimage[2] = xi[2] + dz;
}
/* ----------------------------------------------------------------------
remap the point into the periodic box no matter how far away
adjust 3 image flags encoded in image accordingly
resulting coord must satisfy lo <= coord < hi
MAX is important since coord - prd < lo can happen when coord = hi
for triclinic, point is converted to lamda coords (0-1) before doing remap
image = 10 bits for each dimension
increment/decrement in wrap-around fashion
------------------------------------------------------------------------- */
void Domain::remap(double *x, imageint &image)
{
double *lo,*hi,*period,*coord;
double lamda[3];
imageint idim,otherdims;
if (triclinic == 0) {
lo = boxlo;
hi = boxhi;
period = prd;
coord = x;
} else {
lo = boxlo_lamda;
hi = boxhi_lamda;
period = prd_lamda;
x2lamda(x,lamda);
coord = lamda;
}
if (xperiodic) {
while (coord[0] < lo[0]) {
coord[0] += period[0];
idim = image & IMGMASK;
otherdims = image ^ idim;
idim--;
idim &= IMGMASK;
image = otherdims | idim;
}
while (coord[0] >= hi[0]) {
coord[0] -= period[0];
idim = image & IMGMASK;
otherdims = image ^ idim;
idim++;
idim &= IMGMASK;
image = otherdims | idim;
}
coord[0] = MAX(coord[0],lo[0]);
}
if (yperiodic) {
while (coord[1] < lo[1]) {
coord[1] += period[1];
idim = (image >> IMGBITS) & IMGMASK;
otherdims = image ^ (idim << IMGBITS);
idim--;
idim &= IMGMASK;
image = otherdims | (idim << IMGBITS);
}
while (coord[1] >= hi[1]) {
coord[1] -= period[1];
idim = (image >> IMGBITS) & IMGMASK;
otherdims = image ^ (idim << IMGBITS);
idim++;
idim &= IMGMASK;
image = otherdims | (idim << IMGBITS);
}
coord[1] = MAX(coord[1],lo[1]);
}
if (zperiodic) {
while (coord[2] < lo[2]) {
coord[2] += period[2];
idim = image >> IMG2BITS;
otherdims = image ^ (idim << IMG2BITS);
idim--;
idim &= IMGMASK;
image = otherdims | (idim << IMG2BITS);
}
while (coord[2] >= hi[2]) {
coord[2] -= period[2];
idim = image >> IMG2BITS;
otherdims = image ^ (idim << IMG2BITS);
idim++;
idim &= IMGMASK;
image = otherdims | (idim << IMG2BITS);
}
coord[2] = MAX(coord[2],lo[2]);
}
if (triclinic) lamda2x(coord,x);
}
/* ----------------------------------------------------------------------
remap the point into the periodic box no matter how far away
no image flag calculation
resulting coord must satisfy lo <= coord < hi
MAX is important since coord - prd < lo can happen when coord = hi
for triclinic, point is converted to lamda coords (0-1) before remap
------------------------------------------------------------------------- */
void Domain::remap(double *x)
{
double *lo,*hi,*period,*coord;
double lamda[3];
if (triclinic == 0) {
lo = boxlo;
hi = boxhi;
period = prd;
coord = x;
} else {
lo = boxlo_lamda;
hi = boxhi_lamda;
period = prd_lamda;
x2lamda(x,lamda);
coord = lamda;
}
if (xperiodic) {
while (coord[0] < lo[0]) coord[0] += period[0];
while (coord[0] >= hi[0]) coord[0] -= period[0];
coord[0] = MAX(coord[0],lo[0]);
}
if (yperiodic) {
while (coord[1] < lo[1]) coord[1] += period[1];
while (coord[1] >= hi[1]) coord[1] -= period[1];
coord[1] = MAX(coord[1],lo[1]);
}
if (zperiodic) {
while (coord[2] < lo[2]) coord[2] += period[2];
while (coord[2] >= hi[2]) coord[2] -= period[2];
coord[2] = MAX(coord[2],lo[2]);
}
if (triclinic) lamda2x(coord,x);
}
/* ----------------------------------------------------------------------
remap xnew to be within half box length of xold
do it directly, not iteratively, in case is far away
for triclinic, both points are converted to lamda coords (0-1) before remap
------------------------------------------------------------------------- */
void Domain::remap_near(double *xnew, double *xold)
{
int n;
double *coordnew,*coordold,*period,*half;
double lamdanew[3],lamdaold[3];
if (triclinic == 0) {
period = prd;
half = prd_half;
coordnew = xnew;
coordold = xold;
} else {
period = prd_lamda;
half = prd_half_lamda;
x2lamda(xnew,lamdanew);
coordnew = lamdanew;
x2lamda(xold,lamdaold);
coordold = lamdaold;
}
// iterative form
// if (xperiodic) {
// while (coordnew[0]-coordold[0] > half[0]) coordnew[0] -= period[0];
// while (coordold[0]-coordnew[0] > half[0]) coordnew[0] += period[0];
// }
if (xperiodic) {
if (coordnew[0]-coordold[0] > period[0]) {
n = static_cast<int> ((coordnew[0]-coordold[0])/period[0]);
coordnew[0] -= n*period[0];
}
while (coordnew[0]-coordold[0] > half[0]) coordnew[0] -= period[0];
if (coordold[0]-coordnew[0] > period[0]) {
n = static_cast<int> ((coordold[0]-coordnew[0])/period[0]);
coordnew[0] += n*period[0];
}
while (coordold[0]-coordnew[0] > half[0]) coordnew[0] += period[0];
}
if (yperiodic) {
if (coordnew[1]-coordold[1] > period[1]) {
n = static_cast<int> ((coordnew[1]-coordold[1])/period[1]);
coordnew[1] -= n*period[1];
}
while (coordnew[1]-coordold[1] > half[1]) coordnew[1] -= period[1];
if (coordold[1]-coordnew[1] > period[1]) {
n = static_cast<int> ((coordold[1]-coordnew[1])/period[1]);
coordnew[1] += n*period[1];
}
while (coordold[1]-coordnew[1] > half[1]) coordnew[1] += period[1];
}
if (zperiodic) {
if (coordnew[2]-coordold[2] > period[2]) {
n = static_cast<int> ((coordnew[2]-coordold[2])/period[2]);
coordnew[2] -= n*period[2];
}
while (coordnew[2]-coordold[2] > half[2]) coordnew[2] -= period[2];
if (coordold[2]-coordnew[2] > period[2]) {
n = static_cast<int> ((coordold[2]-coordnew[2])/period[2]);
coordnew[2] += n*period[2];
}
while (coordold[2]-coordnew[2] > half[2]) coordnew[2] += period[2];
}
if (triclinic) lamda2x(coordnew,xnew);
}
/* ----------------------------------------------------------------------
unmap the point via image flags
x overwritten with result, don't reset image flag
for triclinic, use h[] to add in tilt factors in other dims as needed
------------------------------------------------------------------------- */
void Domain::unmap(double *x, imageint image)
{
int xbox = (image & IMGMASK) - IMGMAX;
int ybox = (image >> IMGBITS & IMGMASK) - IMGMAX;
int zbox = (image >> IMG2BITS) - IMGMAX;
if (triclinic == 0) {
x[0] += xbox*xprd;
x[1] += ybox*yprd;
x[2] += zbox*zprd;
} else {
x[0] += h[0]*xbox + h[5]*ybox + h[4]*zbox;
x[1] += h[1]*ybox + h[3]*zbox;
x[2] += h[2]*zbox;
}
}
/* ----------------------------------------------------------------------
unmap the point via image flags
result returned in y, don't reset image flag
for triclinic, use h[] to add in tilt factors in other dims as needed
------------------------------------------------------------------------- */
void Domain::unmap(double *x, imageint image, double *y)
{
int xbox = (image & IMGMASK) - IMGMAX;
int ybox = (image >> IMGBITS & IMGMASK) - IMGMAX;
int zbox = (image >> IMG2BITS) - IMGMAX;
if (triclinic == 0) {
y[0] = x[0] + xbox*xprd;
y[1] = x[1] + ybox*yprd;
y[2] = x[2] + zbox*zprd;
} else {
y[0] = x[0] + h[0]*xbox + h[5]*ybox + h[4]*zbox;
y[1] = x[1] + h[1]*ybox + h[3]*zbox;
y[2] = x[2] + h[2]*zbox;
}
}
/* ----------------------------------------------------------------------
adjust image flags due to triclinic box flip
flip operation is changing box vectors A,B,C to new A',B',C'
A' = A (A does not change)
B' = B + mA (B shifted by A)
C' = C + pB + nA (C shifted by B and/or A)
this requires the image flags change from (a,b,c) to (a',b',c')
so that x_unwrap for each atom is same before/after
x_unwrap_before = xlocal + aA + bB + cC
x_unwrap_after = xlocal + a'A' + b'B' + c'C'
this requires:
c' = c
b' = b - cp
a' = a - (b-cp)m - cn = a - b'm - cn
in other words, for xy flip, change in x flag depends on current y flag
this is b/c the xy flip dramatically changes which tiled image of
simulation box an unwrapped point maps to
------------------------------------------------------------------------- */
void Domain::image_flip(int m, int n, int p)
{
imageint *image = atom->image;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
int xbox = (image[i] & IMGMASK) - IMGMAX;
int ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
int zbox = (image[i] >> IMG2BITS) - IMGMAX;
ybox -= p*zbox;
xbox -= m*ybox + n*zbox;
image[i] = ((imageint) (xbox + IMGMAX) & IMGMASK) |
(((imageint) (ybox + IMGMAX) & IMGMASK) << IMGBITS) |
(((imageint) (zbox + IMGMAX) & IMGMASK) << IMG2BITS);
}
}
/* ----------------------------------------------------------------------
create a lattice
------------------------------------------------------------------------- */
void Domain::set_lattice(int narg, char **arg)
{
if (lattice) delete lattice;
lattice = new Lattice(lmp,narg,arg);
}
/* ----------------------------------------------------------------------
create a new region
------------------------------------------------------------------------- */
void Domain::add_region(int narg, char **arg)
{
if (narg < 2) error->all(FLERR,"Illegal region command");
if (strcmp(arg[1],"delete") == 0) {
delete_region(narg,arg);
return;
}
if (find_region(arg[0]) >= 0) error->all(FLERR,"Reuse of region ID");
// extend Region list if necessary
if (nregion == maxregion) {
maxregion += DELTAREGION;
regions = (Region **)
memory->srealloc(regions,maxregion*sizeof(Region *),"domain:regions");
}
// create the Region
if (lmp->suffix_enable) {
if (lmp->suffix) {
char estyle[256];
sprintf(estyle,"%s/%s",arg[1],lmp->suffix);
-
- if (0) return;
-
-#define REGION_CLASS
-#define RegionStyle(key,Class) \
- else if (strcmp(estyle,#key) == 0) { \
- regions[nregion] = new Class(lmp,narg,arg); \
- regions[nregion]->init(); \
- nregion++; \
- return; \
+ if (region_map->find(estyle) != region_map->end()) {
+ RegionCreator region_creator = (*region_map)[estyle];
+ regions[nregion] = region_creator(lmp, narg, arg);
+ regions[nregion]->init();
+ nregion++;
+ return;
}
-#include "style_region.h"
-#undef RegionStyle
-#undef REGION_CLASS
}
if (lmp->suffix2) {
char estyle[256];
sprintf(estyle,"%s/%s",arg[1],lmp->suffix2);
-
- if (0) return;
-
-#define REGION_CLASS
-#define RegionStyle(key,Class) \
- else if (strcmp(estyle,#key) == 0) { \
- regions[nregion] = new Class(lmp,narg,arg); \
- regions[nregion]->init(); \
- nregion++; \
- return; \
+ if (region_map->find(estyle) != region_map->end()) {
+ RegionCreator region_creator = (*region_map)[estyle];
+ regions[nregion] = region_creator(lmp, narg, arg);
+ regions[nregion]->init();
+ nregion++;
+ return;
}
-#include "style_region.h"
-#undef RegionStyle
-#undef REGION_CLASS
}
}
if (strcmp(arg[1],"none") == 0) error->all(FLERR,"Unknown region style");
-
-#define REGION_CLASS
-#define RegionStyle(key,Class) \
- else if (strcmp(arg[1],#key) == 0) \
- regions[nregion] = new Class(lmp,narg,arg);
-#include "style_region.h"
-#undef REGION_CLASS
-
+ if (region_map->find(arg[1]) != region_map->end()) {
+ RegionCreator region_creator = (*region_map)[arg[1]];
+ regions[nregion] = region_creator(lmp, narg, arg);
+ }
else error->all(FLERR,"Unknown region style");
// initialize any region variables via init()
// in case region is used between runs, e.g. to print a variable
regions[nregion]->init();
nregion++;
}
+/* ----------------------------------------------------------------------
+ one instance per region style in style_region.h
+------------------------------------------------------------------------- */
+
+template <typename T>
+Region *Domain::region_creator(LAMMPS *lmp, int narg, char ** arg)
+{
+ return new T(lmp, narg, arg);
+}
+
/* ----------------------------------------------------------------------
delete a region
------------------------------------------------------------------------- */
void Domain::delete_region(int narg, char **arg)
{
if (narg != 2) error->all(FLERR,"Illegal region command");
int iregion = find_region(arg[0]);
if (iregion == -1) error->all(FLERR,"Delete region ID does not exist");
delete regions[iregion];
regions[iregion] = regions[nregion-1];
nregion--;
}
/* ----------------------------------------------------------------------
return region index if name matches existing region ID
return -1 if no such region
------------------------------------------------------------------------- */
int Domain::find_region(char *name)
{
for (int iregion = 0; iregion < nregion; iregion++)
if (strcmp(name,regions[iregion]->id) == 0) return iregion;
return -1;
}
/* ----------------------------------------------------------------------
(re)set boundary settings
flag = 0, called from the input script
flag = 1, called from change box command
------------------------------------------------------------------------- */
void Domain::set_boundary(int narg, char **arg, int flag)
{
if (narg != 3) error->all(FLERR,"Illegal boundary command");
char c;
for (int idim = 0; idim < 3; idim++)
for (int iside = 0; iside < 2; iside++) {
if (iside == 0) c = arg[idim][0];
else if (iside == 1 && strlen(arg[idim]) == 1) c = arg[idim][0];
else c = arg[idim][1];
if (c == 'p') boundary[idim][iside] = 0;
else if (c == 'f') boundary[idim][iside] = 1;
else if (c == 's') boundary[idim][iside] = 2;
else if (c == 'm') boundary[idim][iside] = 3;
else {
if (flag == 0) error->all(FLERR,"Illegal boundary command");
if (flag == 1) error->all(FLERR,"Illegal change_box command");
}
}
for (int idim = 0; idim < 3; idim++)
if ((boundary[idim][0] == 0 && boundary[idim][1]) ||
(boundary[idim][0] && boundary[idim][1] == 0))
error->all(FLERR,"Both sides of boundary must be periodic");
if (boundary[0][0] == 0) xperiodic = 1;
else xperiodic = 0;
if (boundary[1][0] == 0) yperiodic = 1;
else yperiodic = 0;
if (boundary[2][0] == 0) zperiodic = 1;
else zperiodic = 0;
periodicity[0] = xperiodic;
periodicity[1] = yperiodic;
periodicity[2] = zperiodic;
nonperiodic = 0;
if (xperiodic == 0 || yperiodic == 0 || zperiodic == 0) {
nonperiodic = 1;
if (boundary[0][0] >= 2 || boundary[0][1] >= 2 ||
boundary[1][0] >= 2 || boundary[1][1] >= 2 ||
boundary[2][0] >= 2 || boundary[2][1] >= 2) nonperiodic = 2;
}
}
/* ----------------------------------------------------------------------
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
------------------------------------------------------------------------- */
void Domain::print_box(const char *str)
{
if (comm->me == 0) {
if (screen) {
if (triclinic == 0)
fprintf(screen,"%sorthogonal box = (%g %g %g) to (%g %g %g)\n",
str,boxlo[0],boxlo[1],boxlo[2],boxhi[0],boxhi[1],boxhi[2]);
else {
char *format = (char *)
"%striclinic box = (%g %g %g) to (%g %g %g) with tilt (%g %g %g)\n";
fprintf(screen,format,
str,boxlo[0],boxlo[1],boxlo[2],boxhi[0],boxhi[1],boxhi[2],
xy,xz,yz);
}
}
if (logfile) {
if (triclinic == 0)
fprintf(logfile,"%sorthogonal box = (%g %g %g) to (%g %g %g)\n",
str,boxlo[0],boxlo[1],boxlo[2],boxhi[0],boxhi[1],boxhi[2]);
else {
char *format = (char *)
"%striclinic box = (%g %g %g) to (%g %g %g) with tilt (%g %g %g)\n";
fprintf(logfile,format,
str,boxlo[0],boxlo[1],boxlo[2],boxhi[0],boxhi[1],boxhi[2],
xy,xz,yz);
}
}
}
}
/* ----------------------------------------------------------------------
format boundary string for output
assume str is 9 chars or more in length
------------------------------------------------------------------------- */
void Domain::boundary_string(char *str)
{
int m = 0;
for (int idim = 0; idim < 3; idim++) {
for (int iside = 0; iside < 2; iside++) {
if (boundary[idim][iside] == 0) str[m++] = 'p';
else if (boundary[idim][iside] == 1) str[m++] = 'f';
else if (boundary[idim][iside] == 2) str[m++] = 's';
else if (boundary[idim][iside] == 3) str[m++] = 'm';
}
str[m++] = ' ';
}
str[8] = '\0';
}
/* ----------------------------------------------------------------------
convert triclinic 0-1 lamda coords to box coords for all N atoms
x = H lamda + x0;
------------------------------------------------------------------------- */
void Domain::lamda2x(int n)
{
double **x = atom->x;
for (int i = 0; i < n; i++) {
x[i][0] = h[0]*x[i][0] + h[5]*x[i][1] + h[4]*x[i][2] + boxlo[0];
x[i][1] = h[1]*x[i][1] + h[3]*x[i][2] + boxlo[1];
x[i][2] = h[2]*x[i][2] + boxlo[2];
}
}
/* ----------------------------------------------------------------------
convert box coords to triclinic 0-1 lamda coords for all N atoms
lamda = H^-1 (x - x0)
------------------------------------------------------------------------- */
void Domain::x2lamda(int n)
{
double delta[3];
double **x = atom->x;
for (int i = 0; i < n; i++) {
delta[0] = x[i][0] - boxlo[0];
delta[1] = x[i][1] - boxlo[1];
delta[2] = x[i][2] - boxlo[2];
x[i][0] = h_inv[0]*delta[0] + h_inv[5]*delta[1] + h_inv[4]*delta[2];
x[i][1] = h_inv[1]*delta[1] + h_inv[3]*delta[2];
x[i][2] = h_inv[2]*delta[2];
}
}
/* ----------------------------------------------------------------------
convert triclinic 0-1 lamda coords to box coords for one atom
x = H lamda + x0;
lamda and x can point to same 3-vector
------------------------------------------------------------------------- */
void Domain::lamda2x(double *lamda, double *x)
{
x[0] = h[0]*lamda[0] + h[5]*lamda[1] + h[4]*lamda[2] + boxlo[0];
x[1] = h[1]*lamda[1] + h[3]*lamda[2] + boxlo[1];
x[2] = h[2]*lamda[2] + boxlo[2];
}
/* ----------------------------------------------------------------------
convert box coords to triclinic 0-1 lamda coords for one atom
lamda = H^-1 (x - x0)
x and lamda can point to same 3-vector
------------------------------------------------------------------------- */
void Domain::x2lamda(double *x, double *lamda)
{
double delta[3];
delta[0] = x[0] - boxlo[0];
delta[1] = x[1] - boxlo[1];
delta[2] = x[2] - boxlo[2];
lamda[0] = h_inv[0]*delta[0] + h_inv[5]*delta[1] + h_inv[4]*delta[2];
lamda[1] = h_inv[1]*delta[1] + h_inv[3]*delta[2];
lamda[2] = h_inv[2]*delta[2];
}
/* ----------------------------------------------------------------------
convert box coords to triclinic 0-1 lamda coords for one atom
use my_boxlo & my_h_inv stored by caller for previous state of box
lamda = H^-1 (x - x0)
x and lamda can point to same 3-vector
------------------------------------------------------------------------- */
void Domain::x2lamda(double *x, double *lamda,
double *my_boxlo, double *my_h_inv)
{
double delta[3];
delta[0] = x[0] - my_boxlo[0];
delta[1] = x[1] - my_boxlo[1];
delta[2] = x[2] - my_boxlo[2];
lamda[0] = my_h_inv[0]*delta[0] + my_h_inv[5]*delta[1] + my_h_inv[4]*delta[2];
lamda[1] = my_h_inv[1]*delta[1] + my_h_inv[3]*delta[2];
lamda[2] = my_h_inv[2]*delta[2];
}
/* ----------------------------------------------------------------------
convert 8 lamda corner pts of lo/hi box to box coords
return bboxlo/hi = bounding box around 8 corner pts in box coords
------------------------------------------------------------------------- */
void Domain::bbox(double *lo, double *hi, double *bboxlo, double *bboxhi)
{
double x[3];
bboxlo[0] = bboxlo[1] = bboxlo[2] = BIG;
bboxhi[0] = bboxhi[1] = bboxhi[2] = -BIG;
x[0] = lo[0]; x[1] = lo[1]; x[2] = lo[2];
lamda2x(x,x);
bboxlo[0] = MIN(bboxlo[0],x[0]); bboxhi[0] = MAX(bboxhi[0],x[0]);
bboxlo[1] = MIN(bboxlo[1],x[1]); bboxhi[1] = MAX(bboxhi[1],x[1]);
bboxlo[2] = MIN(bboxlo[2],x[2]); bboxhi[2] = MAX(bboxhi[2],x[2]);
x[0] = hi[0]; x[1] = lo[1]; x[2] = lo[2];
lamda2x(x,x);
bboxlo[0] = MIN(bboxlo[0],x[0]); bboxhi[0] = MAX(bboxhi[0],x[0]);
bboxlo[1] = MIN(bboxlo[1],x[1]); bboxhi[1] = MAX(bboxhi[1],x[1]);
bboxlo[2] = MIN(bboxlo[2],x[2]); bboxhi[2] = MAX(bboxhi[2],x[2]);
x[0] = lo[0]; x[1] = hi[1]; x[2] = lo[2];
lamda2x(x,x);
bboxlo[0] = MIN(bboxlo[0],x[0]); bboxhi[0] = MAX(bboxhi[0],x[0]);
bboxlo[1] = MIN(bboxlo[1],x[1]); bboxhi[1] = MAX(bboxhi[1],x[1]);
bboxlo[2] = MIN(bboxlo[2],x[2]); bboxhi[2] = MAX(bboxhi[2],x[2]);
x[0] = hi[0]; x[1] = hi[1]; x[2] = lo[2];
lamda2x(x,x);
bboxlo[0] = MIN(bboxlo[0],x[0]); bboxhi[0] = MAX(bboxhi[0],x[0]);
bboxlo[1] = MIN(bboxlo[1],x[1]); bboxhi[1] = MAX(bboxhi[1],x[1]);
bboxlo[2] = MIN(bboxlo[2],x[2]); bboxhi[2] = MAX(bboxhi[2],x[2]);
x[0] = lo[0]; x[1] = lo[1]; x[2] = hi[2];
lamda2x(x,x);
bboxlo[0] = MIN(bboxlo[0],x[0]); bboxhi[0] = MAX(bboxhi[0],x[0]);
bboxlo[1] = MIN(bboxlo[1],x[1]); bboxhi[1] = MAX(bboxhi[1],x[1]);
bboxlo[2] = MIN(bboxlo[2],x[2]); bboxhi[2] = MAX(bboxhi[2],x[2]);
x[0] = hi[0]; x[1] = lo[1]; x[2] = hi[2];
lamda2x(x,x);
bboxlo[0] = MIN(bboxlo[0],x[0]); bboxhi[0] = MAX(bboxhi[0],x[0]);
bboxlo[1] = MIN(bboxlo[1],x[1]); bboxhi[1] = MAX(bboxhi[1],x[1]);
bboxlo[2] = MIN(bboxlo[2],x[2]); bboxhi[2] = MAX(bboxhi[2],x[2]);
x[0] = lo[0]; x[1] = hi[1]; x[2] = hi[2];
lamda2x(x,x);
bboxlo[0] = MIN(bboxlo[0],x[0]); bboxhi[0] = MAX(bboxhi[0],x[0]);
bboxlo[1] = MIN(bboxlo[1],x[1]); bboxhi[1] = MAX(bboxhi[1],x[1]);
bboxlo[2] = MIN(bboxlo[2],x[2]); bboxhi[2] = MAX(bboxhi[2],x[2]);
x[0] = hi[0]; x[1] = hi[1]; x[2] = hi[2];
lamda2x(x,x);
bboxlo[0] = MIN(bboxlo[0],x[0]); bboxhi[0] = MAX(bboxhi[0],x[0]);
bboxlo[1] = MIN(bboxlo[1],x[1]); bboxhi[1] = MAX(bboxhi[1],x[1]);
bboxlo[2] = MIN(bboxlo[2],x[2]); bboxhi[2] = MAX(bboxhi[2],x[2]);
}
/* ----------------------------------------------------------------------
compute 8 corner pts of my triclinic sub-box
output is in corners, see ordering in lamda_box_corners
------------------------------------------------------------------------- */
void Domain::box_corners()
{
lamda_box_corners(boxlo_lamda,boxhi_lamda);
}
/* ----------------------------------------------------------------------
compute 8 corner pts of my triclinic sub-box
output is in corners, see ordering in lamda_box_corners
------------------------------------------------------------------------- */
void Domain::subbox_corners()
{
lamda_box_corners(sublo_lamda,subhi_lamda);
}
/* ----------------------------------------------------------------------
compute 8 corner pts of any triclinic box with lo/hi in lamda coords
8 output conners are ordered with x changing fastest, then y, finally z
could be more efficient if just coded with xy,yz,xz explicitly
------------------------------------------------------------------------- */
void Domain::lamda_box_corners(double *lo, double *hi)
{
corners[0][0] = lo[0]; corners[0][1] = lo[1]; corners[0][2] = lo[2];
lamda2x(corners[0],corners[0]);
corners[1][0] = hi[0]; corners[1][1] = lo[1]; corners[1][2] = lo[2];
lamda2x(corners[1],corners[1]);
corners[2][0] = lo[0]; corners[2][1] = hi[1]; corners[2][2] = lo[2];
lamda2x(corners[2],corners[2]);
corners[3][0] = hi[0]; corners[3][1] = hi[1]; corners[3][2] = lo[2];
lamda2x(corners[3],corners[3]);
corners[4][0] = lo[0]; corners[4][1] = lo[1]; corners[4][2] = hi[2];
lamda2x(corners[4],corners[4]);
corners[5][0] = hi[0]; corners[5][1] = lo[1]; corners[5][2] = hi[2];
lamda2x(corners[5],corners[5]);
corners[6][0] = lo[0]; corners[6][1] = hi[1]; corners[6][2] = hi[2];
lamda2x(corners[6],corners[6]);
corners[7][0] = hi[0]; corners[7][1] = hi[1]; corners[7][2] = subhi_lamda[2];
lamda2x(corners[7],corners[7]);
}
diff --git a/src/domain.h b/src/domain.h
index 8654b5916..ef3de4207 100644
--- a/src/domain.h
+++ b/src/domain.h
@@ -1,269 +1,278 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifndef LMP_DOMAIN_H
#define LMP_DOMAIN_H
#include <math.h>
#include "pointers.h"
+#include <map>
+#include <string>
namespace LAMMPS_NS {
class Domain : protected Pointers {
public:
int box_exist; // 0 = not yet created, 1 = exists
int dimension; // 2 = 2d, 3 = 3d
int nonperiodic; // 0 = periodic in all 3 dims
// 1 = periodic or fixed in all 6
// 2 = shrink-wrap in any of 6
int xperiodic,yperiodic,zperiodic; // 0 = non-periodic, 1 = periodic
int periodicity[3]; // xyz periodicity as array
int boundary[3][2]; // settings for 6 boundaries
// 0 = periodic
// 1 = fixed non-periodic
// 2 = shrink-wrap non-periodic
// 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
double xprd_half,yprd_half,zprd_half; // half dimensions
double prd[3]; // array form of dimensions
double prd_half[3]; // array form of half dimensions
// triclinic box
// xprd,xprd_half,prd,prd_half =
// same as if untilted
double prd_lamda[3]; // lamda box = (1,1,1)
double prd_half_lamda[3]; // lamda half box = (0.5,0.5,0.5)
double boxlo[3],boxhi[3]; // orthogonal box global bounds
// triclinic box
// boxlo/hi = same as if untilted
double boxlo_lamda[3],boxhi_lamda[3]; // lamda box = (0,1)
double boxlo_bound[3],boxhi_bound[3]; // bounding box of tilted domain
double corners[8][3]; // 8 corner points
// orthogonal box & triclinic box
double minxlo,minxhi; // minimum size of global box
double minylo,minyhi; // when shrink-wrapping
double minzlo,minzhi; // tri only possible for non-skew dims
// orthogonal box
double sublo[3],subhi[3]; // sub-box bounds on this proc
// triclinic box
// sublo/hi = undefined
double sublo_lamda[3],subhi_lamda[3]; // bounds of subbox in lamda
// triclinic box
double xy,xz,yz; // 3 tilt factors
double h[6],h_inv[6]; // shape matrix in Voigt notation
double h_rate[6],h_ratelo[3]; // rate of box size/shape change
int box_change; // 1 if any of next 3 flags are set, else 0
int box_change_size; // 1 if box size changes, 0 if not
int box_change_shape; // 1 if box shape changes, 0 if not
int box_change_domain; // 1 if proc sub-domains change, 0 if not
int deform_flag; // 1 if fix deform exist, else 0
int deform_vremap; // 1 if fix deform remaps v, else 0
int deform_groupbit; // atom group to perform v remap for
class Lattice *lattice; // user-defined lattice
int nregion; // # of defined Regions
int maxregion; // max # list can hold
class Region **regions; // list of defined Regions
int copymode;
+ typedef Region *(*RegionCreator)(LAMMPS *,int,char**);
+ typedef std::map<std::string,RegionCreator> RegionCreatorMap;
+ RegionCreatorMap *region_map;
+
Domain(class LAMMPS *);
virtual ~Domain();
virtual void init();
void set_initial_box(int expandflag=1);
virtual void set_global_box();
virtual void set_lamda_box();
virtual void set_local_box();
virtual void reset_box();
virtual void pbc();
void image_check();
void box_too_small_check();
void subbox_too_small_check(double);
void minimum_image(double &, double &, double &);
void minimum_image(double *);
int closest_image(int, int);
void closest_image(const double * const, const double * const,
double * const);
void remap(double *, imageint &);
void remap(double *);
void remap_near(double *, double *);
void unmap(double *, imageint);
void unmap(double *, imageint, double *);
void image_flip(int, int, int);
void set_lattice(int, char **);
void add_region(int, char **);
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 *);
virtual void lamda2x(int);
virtual void x2lamda(int);
virtual void lamda2x(double *, double *);
virtual void x2lamda(double *, double *);
int inside(double *);
int inside_nonperiodic(double *);
void x2lamda(double *, double *, double *, double *);
void bbox(double *, double *, double *, double *);
void box_corners();
void subbox_corners();
void lamda_box_corners(double *, double *);
// minimum image convention check
// return 1 if any distance > 1/2 of box size
// indicates a special neighbor is actually not in a bond,
// but is a far-away image that should be treated as an unbonded neighbor
// inline since called from neighbor build inner loop
//
inline int minimum_image_check(double dx, double dy, double dz) {
if (xperiodic && fabs(dx) > xprd_half) return 1;
if (yperiodic && fabs(dy) > yprd_half) return 1;
if (zperiodic && fabs(dz) > zprd_half) return 1;
return 0;
}
protected:
double small[3]; // fractions of box lengths
+
+ private:
+ template <typename T> static Region *region_creator(LAMMPS *,int,char**);
};
}
#endif
/* ERROR/WARNING messages:
E: Box bounds are invalid
The box boundaries specified in the read_data file are invalid. The
lo value must be less than the hi value for all 3 dimensions.
E: Cannot skew triclinic box in z for 2d simulation
Self-explanatory.
E: Triclinic box skew is too large
The displacement in a skewed direction must be less than half the box
length in that dimension. E.g. the xy tilt must be between -half and
+half of the x box length. This constraint can be relaxed by using
the box tilt command.
W: Triclinic box skew is large
The displacement in a skewed direction is normally required to be less
than half the box length in that dimension. E.g. the xy tilt must be
between -half and +half of the x box length. You have relaxed the
constraint using the box tilt command, but the warning means that a
LAMMPS simulation may be inefficient as a result.
E: Illegal simulation box
The lower bound of the simulation box is greater than the upper bound.
E: Bond atom missing in image check
The 2nd atom in a particular bond is missing on this processor.
Typically this is because the pairwise cutoff is set too short or the
bond has blown apart and an atom is too far away.
W: Inconsistent image flags
The image flags for a pair on bonded atoms appear to be inconsistent.
Inconsistent means that when the coordinates of the two atoms are
unwrapped using the image flags, the two atoms are far apart.
Specifically they are further apart than half a periodic box length.
Or they are more than a box length apart in a non-periodic dimension.
This is usually due to the initial data file not having correct image
flags for the 2 atoms in a bond that straddles a periodic boundary.
They should be different by 1 in that case. This is a warning because
inconsistent image flags will not cause problems for dynamics or most
LAMMPS simulations. However they can cause problems when such atoms
are used with the fix rigid or replicate commands.
W: Bond atom missing in image check
The 2nd atom in a particular bond is missing on this processor.
Typically this is because the pairwise cutoff is set too short or the
bond has blown apart and an atom is too far away.
E: Bond atom missing in box size check
The 2nd atoms needed to compute a particular bond is missing on this
processor. Typically this is because the pairwise cutoff is set too
short or the bond has blown apart and an atom is too far away.
W: Bond atom missing in box size check
The 2nd atoms needed to compute a particular bond is missing on this
processor. Typically this is because the pairwise cutoff is set too
short or the bond has blown apart and an atom is too far away.
W: Bond/angle/dihedral extent > half of periodic box length
This is a restriction because LAMMPS can be confused about which image
of an atom in the bonded interaction is the correct one to use.
"Extent" in this context means the maximum end-to-end length of the
bond/angle/dihedral. LAMMPS computes this by taking the maximum bond
length, multiplying by the number of bonds in the interaction (e.g. 3
for a dihedral) and adding a small amount of stretch.
W: Proc sub-domain size < neighbor skin, could lead to lost atoms
The decomposition of the physical domain (likely due to load
balancing) has led to a processor's sub-domain being smaller than the
neighbor skin in one or more dimensions. Since reneighboring is
triggered by atoms moving the skin distance, this may lead to lost
atoms, if an atom moves all the way across a neighboring processor's
sub-domain before reneighboring is triggered.
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Reuse of region ID
A region ID cannot be used twice.
E: Unknown region style
The choice of region style is unknown.
E: Delete region ID does not exist
Self-explanatory.
E: Both sides of boundary must be periodic
Cannot specify a boundary as periodic only on the lo or hi side. Must
be periodic on both sides.
*/
diff --git a/src/info.cpp b/src/info.cpp
index 22cb432ac..c150e0a88 100644
--- a/src/info.cpp
+++ b/src/info.cpp
@@ -1,1018 +1,1019 @@
/* ----------------------------------------------------------------------
fputs
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors: Axel Kohlmeyer (Temple U),
Richard Berger (Temple U)
------------------------------------------------------------------------- */
#include <string.h>
#include "info.h"
#include "accelerator_kokkos.h"
#include "atom.h"
#include "comm.h"
#include "compute.h"
#include "domain.h"
#include "dump.h"
#include "fix.h"
#include "force.h"
#include "pair.h"
#include "pair_hybrid.h"
#include "group.h"
#include "input.h"
#include "modify.h"
#include "neighbor.h"
#include "output.h"
#include "region.h"
#include "universe.h"
#include "variable.h"
#include "update.h"
#include "error.h"
#include <time.h>
#include <vector>
#include <string>
#include <algorithm>
#ifdef _WIN32
#define PSAPI_VERSION=1
#include <windows.h>
#include <stdint.h>
#include <psapi.h>
#else
#include <sys/time.h>
#include <sys/resource.h>
#include <sys/utsname.h>
#endif
#if defined __linux
#include <malloc.h>
#endif
namespace LAMMPS_NS {
// same as in variable.cpp
enum {INDEX,LOOP,WORLD,UNIVERSE,ULOOP,STRING,GETENV,
SCALARFILE,ATOMFILE,FORMAT,EQUAL,ATOM,PYTHON};
enum {COMPUTES=1<<0,
DUMPS=1<<1,
FIXES=1<<2,
GROUPS=1<<3,
REGIONS=1<<4,
CONFIG=1<<5,
TIME=1<<6,
VARIABLES=1<<7,
SYSTEM=1<<8,
COMM=1<<9,
ATOM_STYLES=1<<10,
INTEGRATE_STYLES=1<<11,
MINIMIZE_STYLES=1<<12,
PAIR_STYLES=1<<13,
BOND_STYLES=1<<14,
ANGLE_STYLES=1<<15,
DIHEDRAL_STYLES=1<<16,
IMPROPER_STYLES=1<<17,
KSPACE_STYLES=1<<18,
FIX_STYLES=1<<19,
COMPUTE_STYLES=1<<20,
REGION_STYLES=1<<21,
DUMP_STYLES=1<<22,
COMMAND_STYLES=1<<23,
ALL=~0};
static const int STYLES = ATOM_STYLES | INTEGRATE_STYLES | MINIMIZE_STYLES | PAIR_STYLES | BOND_STYLES | \
ANGLE_STYLES | DIHEDRAL_STYLES | IMPROPER_STYLES | KSPACE_STYLES | FIX_STYLES | \
COMPUTE_STYLES | REGION_STYLES | DUMP_STYLES | COMMAND_STYLES;
}
static const char *varstyles[] = {
"index", "loop", "world", "universe", "uloop", "string", "getenv",
"file", "atomfile", "format", "equal", "atom", "python", "(unknown)"};
static const char *mapstyles[] = { "none", "array", "hash" };
static const char *commstyles[] = { "brick", "tiled" };
static const char *commlayout[] = { "uniform", "nonuniform", "irregular" };
static const char bstyles[] = "pfsm";
using namespace LAMMPS_NS;
using namespace std;
static void print_columns(FILE* fp, vector<string> & styles);
/* ---------------------------------------------------------------------- */
void Info::command(int narg, char **arg)
{
FILE *out=screen;
int flags=0;
if (comm->me != 0) return;
// parse arguments
int idx = 0;
while (idx < narg) {
if (strncmp(arg[idx],"all",3) == 0) {
flags |= ALL;
++idx;
} else if ((idx+1 < narg) && (strncmp(arg[idx],"out",3) == 0)
&& (strncmp(arg[idx+1],"screen",3) == 0)) {
if ((out != screen) && (out != logfile)) fclose(out);
out = screen;
idx += 2;
} else if ((idx+1 < narg) && (strncmp(arg[idx],"out",3) == 0)
&& (strncmp(arg[idx+1],"log",3) == 0)) {
if ((out != screen) && (out != logfile)) fclose(out);
out = logfile;
idx += 2;
} else if ((idx+2 < narg) && (strncmp(arg[idx],"out",3) == 0)
&& (strncmp(arg[idx+1],"append",3) == 0)) {
if ((out != screen) && (out != logfile)) fclose(out);
out = fopen(arg[idx+2],"a");
idx += 3;
} else if ((idx+2 < narg) && (strncmp(arg[idx],"out",3) == 0)
&& (strncmp(arg[idx+1],"overwrite",3) == 0)) {
if ((out != screen) && (out != logfile)) fclose(out);
out = fopen(arg[idx+2],"w");
idx += 3;
} else if (strncmp(arg[idx],"atom_styles",3) == 0) {
flags |= ATOM_STYLES;
++idx;
} else if (strncmp(arg[idx],"integrate_styles",3) == 0) {
flags |= INTEGRATE_STYLES;
++idx;
} else if (strncmp(arg[idx],"minimize_styles",3) == 0) {
flags |= MINIMIZE_STYLES;
++idx;
} else if (strncmp(arg[idx],"pair_styles",3) == 0) {
flags |= PAIR_STYLES;
++idx;
} else if (strncmp(arg[idx],"bond_styles",3) == 0) {
flags |= BOND_STYLES;
++idx;
} else if (strncmp(arg[idx],"angle_styles",3) == 0) {
flags |= ANGLE_STYLES;
++idx;
} else if (strncmp(arg[idx],"dihedral_styles",3) == 0) {
flags |= DIHEDRAL_STYLES;
++idx;
} else if (strncmp(arg[idx],"improper_styles",3) == 0) {
flags |= IMPROPER_STYLES;
++idx;
} else if (strncmp(arg[idx],"kspace_styles",3) == 0) {
flags |= KSPACE_STYLES;
++idx;
} else if (strncmp(arg[idx],"fix_styles",3) == 0) {
flags |= FIX_STYLES;
++idx;
} else if (strncmp(arg[idx],"compute_styles",5) == 0) {
flags |= COMPUTE_STYLES;
++idx;
} else if (strncmp(arg[idx],"region_styles",3) == 0) {
flags |= REGION_STYLES;
++idx;
} else if (strncmp(arg[idx],"dump_styles",3) == 0) {
flags |= DUMP_STYLES;
++idx;
} else if (strncmp(arg[idx],"command_styles",5) == 0) {
flags |= COMMAND_STYLES;
++idx;
} else if (strncmp(arg[idx],"communication",5) == 0) {
flags |= COMM;
++idx;
} else if (strncmp(arg[idx],"computes",5) == 0) {
flags |= COMPUTES;
++idx;
} else if (strncmp(arg[idx],"dumps",5) == 0) {
flags |= DUMPS;
++idx;
} else if (strncmp(arg[idx],"fixes",5) == 0) {
flags |= FIXES;
++idx;
} else if (strncmp(arg[idx],"groups",3) == 0) {
flags |= GROUPS;
++idx;
} else if (strncmp(arg[idx],"regions",3) == 0) {
flags |= REGIONS;
++idx;
} else if (strncmp(arg[idx],"config",3) == 0) {
flags |= CONFIG;
++idx;
} else if (strncmp(arg[idx],"time",3) == 0) {
flags |= TIME;
++idx;
} else if (strncmp(arg[idx],"variables",3) == 0) {
flags |= VARIABLES;
++idx;
} else if (strncmp(arg[idx],"system",3) == 0) {
flags |= SYSTEM;
++idx;
} else if (strncmp(arg[idx],"styles",3) == 0) {
flags |= STYLES;
++idx;
} else {
error->warning(FLERR,"Ignoring unknown or incorrect info command flag");
++idx;
}
}
if (out == NULL) return;
fputs("\nInfo-Info-Info-Info-Info-Info-Info-Info-Info-Info-Info\n",out);
time_t now = time(NULL);
fprintf(out,"Printed on %s\n",ctime(&now));
if (flags & CONFIG) {
fprintf(out,"\nLAMMPS version: %s / %s\n",
universe->version, universe->num_ver);
fprintf(out,"sizeof(smallint): %3d-bit\n",(int)sizeof(smallint)*8);
fprintf(out,"sizeof(imageint): %3d-bit\n",(int)sizeof(imageint)*8);
fprintf(out,"sizeof(tagint): %3d-bit\n",(int)sizeof(tagint)*8);
fprintf(out,"sizeof(bigint): %3d-bit\n",(int)sizeof(bigint)*8);
#if defined(_WIN32)
DWORD fullversion,majorv,minorv,buildv=0;
fullversion = GetVersion();
majorv = (DWORD) (LOBYTE(LOWORD(fullversion)));
minorv = (DWORD) (HIBYTE(LOWORD(fullversion)));
if (fullversion < 0x80000000)
buildv = (DWORD) (HIWORD(fullversion));
SYSTEM_INFO si;
GetSystemInfo(&si);
const char *machine;
switch (si.wProcessorArchitecture) {
case PROCESSOR_ARCHITECTURE_AMD64:
machine = (const char *) "x86_64";
break;
case PROCESSOR_ARCHITECTURE_ARM:
machine = (const char *) "arm";
break;
case PROCESSOR_ARCHITECTURE_IA64:
machine = (const char *) "ia64";
break;
case PROCESSOR_ARCHITECTURE_INTEL:
machine = (const char *) "i386";
break;
default:
machine = (const char *) "(unknown)";
}
fprintf(out,"\nOS information: Windows %d.%d (%d) on %s\n",
majorv,minorv,buildv,machine);
#else
struct utsname ut;
uname(&ut);
fprintf(out,"\nOS information: %s %s on %s\n",
ut.sysname, ut.release, ut.machine);
#endif
fprintf(out,"\nMemory allocation information (MPI rank 0)\n");
#if defined(_WIN32)
HANDLE phandle = GetCurrentProcess();
PROCESS_MEMORY_COUNTERS_EX pmc;
GetProcessMemoryInfo(phandle,(PROCESS_MEMORY_COUNTERS *)&pmc,sizeof(pmc));
fprintf(out,"Non-shared memory use: %.3g Mbyte\n",
(double)pmc.PrivateUsage/1048576.0);
fprintf(out,"Maximum working set size: %.3g Mbyte\n",
(double)pmc.PeakWorkingSetSize/1048576.0);
#else
#if defined(__linux)
struct mallinfo mi;
mi = mallinfo();
fprintf(out,"Total dynamically allocated memory: %.3g Mbyte\n",
(double)mi.uordblks/1048576.0);
#endif
struct rusage ru;
if (getrusage(RUSAGE_SELF, &ru) == 0) {
fprintf(out,"Maximum resident set size: %.3g Mbyte\n",
(double)ru.ru_maxrss/1024.0);
}
#endif
}
if (flags & COMM) {
int major,minor;
MPI_Get_version(&major,&minor);
fprintf(out,"\nCommunication information:\n");
fprintf(out,"MPI library level: MPI v%d.%d\n",major,minor);
fprintf(out,"Comm style = %s, Comm layout = %s\n",
commstyles[comm->style], commlayout[comm->layout]);
fprintf(out,"Communicate velocities for ghost atoms = %s\n",
comm->ghost_velocity ? "yes" : "no");
if (comm->mode == 0) {
fprintf(out,"Communication mode = single\n");
fprintf(out,"Communication cutoff = %g\n",
MAX(comm->cutghostuser,neighbor->cutneighmax));
}
if (comm->mode == 1) {
fprintf(out,"Communication mode = multi\n");
double cut;
for (int i=1; i <= atom->ntypes && neighbor->cuttype; ++i) {
cut = neighbor->cuttype[i];
if (comm->cutusermulti) cut = MAX(cut,comm->cutusermulti[i]);
fprintf(out,"Communication cutoff for type %d = %g\n", i, cut);
}
}
fprintf(out,"Nprocs = %d, Nthreads = %d\n",
comm->nprocs, comm->nthreads);
if (domain->box_exist)
fprintf(out,"Processor grid = %d x %d x %d\n",comm->procgrid[0],
comm->procgrid[1], comm->procgrid[2]);
}
if (flags & SYSTEM) {
fprintf(out,"\nSystem information:\n");
fprintf(out,"Units = %s\n",update->unit_style);
fprintf(out,"Atom style = %s\n", atom->atom_style);
fprintf(out,"Atom map = %s\n", mapstyles[atom->map_style]);
if (atom->molecular > 0) {
const char *msg;
msg = (atom->molecular == 2) ? "template" : "standard";
fprintf(out,"Molecule type = %s\n",msg);
}
fprintf(out,"Atoms = " BIGINT_FORMAT ", types = %d, style = %s\n",
atom->natoms, atom->ntypes, force->pair_style);
if (force->pair && strstr(force->pair_style,"hybrid")) {
PairHybrid *hybrid = (PairHybrid *)force->pair;
fprintf(out,"Hybrid sub-styles:");
for (int i=0; i < hybrid->nstyles; ++i)
fprintf(out," %s", hybrid->keywords[i]);
fputc('\n',out);
}
if (atom->molecular > 0) {
const char *msg;
msg = force->bond_style ? force->bond_style : "none";
fprintf(out,"Bonds = " BIGINT_FORMAT ", types = %d, style = %s\n",
atom->nbonds, atom->nbondtypes, msg);
msg = force->angle_style ? force->angle_style : "none";
fprintf(out,"Angles = " BIGINT_FORMAT ", types = %d, style = %s\n",
atom->nangles, atom->nangletypes, msg);
msg = force->dihedral_style ? force->dihedral_style : "none";
fprintf(out,"Dihedrals = " BIGINT_FORMAT ", types = %d, style = %s\n",
atom->ndihedrals, atom->ndihedraltypes, msg);
msg = force->improper_style ? force->improper_style : "none";
fprintf(out,"Impropers = " BIGINT_FORMAT ", types = %d, style = %s\n",
atom->nimpropers, atom->nimpropertypes, msg);
const double * const special_lj = force->special_lj;
const double * const special_coul = force->special_coul;
fprintf(out,"Special bond factors lj = %-10g %-10g %-10g\n"
"Special bond factors coul = %-10g %-10g %-10g\n",
special_lj[1],special_lj[2],special_lj[3],
special_coul[1],special_coul[2],special_coul[3]);
}
fprintf(out,"Kspace style = %s\n",
force->kspace ? force->kspace_style : "none");
if (domain->box_exist) {
fprintf(out,"\nDimensions = %d\n",domain->dimension);
fprintf(out,"%s box = %g x %g x %g\n",
domain->triclinic ? "Triclinic" : "Orthogonal",
domain->xprd, domain->yprd, domain->zprd);
fprintf(out,"Boundaries = %c,%c %c,%c %c,%c\n",
bstyles[domain->boundary[0][0]],bstyles[domain->boundary[0][1]],
bstyles[domain->boundary[1][0]],bstyles[domain->boundary[1][1]],
bstyles[domain->boundary[2][0]],bstyles[domain->boundary[2][1]]);
fprintf(out,"xlo, xhi = %g, %g\n", domain->boxlo[0], domain->boxhi[0]);
fprintf(out,"ylo, yhi = %g, %g\n", domain->boxlo[1], domain->boxhi[1]);
fprintf(out,"zlo, zhi = %g, %g\n", domain->boxlo[2], domain->boxhi[2]);
if (domain->triclinic)
fprintf(out,"Xy, xz, yz = %g, %g, %g\n",
domain->xy, domain->xz, domain->yz);
} else {
fputs("\nBox has not yet been created\n",out);
}
}
if (flags & GROUPS) {
int ngroup = group->ngroup;
char **names = group->names;
int *dynamic = group->dynamic;
fprintf(out,"\nGroup information:\n");
for (int i=0; i < ngroup; ++i) {
if (names[i])
fprintf(out,"Group[%2d]: %s (%s)\n",
i, names[i], dynamic[i] ? "dynamic" : "static");
}
}
if (flags & REGIONS) {
int nreg = domain->nregion;
Region **regs = domain->regions;
fprintf(out,"\nRegion information:\n");
for (int i=0; i < nreg; ++i) {
fprintf(out,"Region[%3d]: %s, style = %s, side = %s\n",
i, regs[i]->id, regs[i]->style,
regs[i]->interior ? "in" : "out");
}
}
if (flags & COMPUTES) {
int ncompute = modify->ncompute;
Compute **compute = modify->compute;
char **names = group->names;
fprintf(out,"\nCompute information:\n");
for (int i=0; i < ncompute; ++i) {
fprintf(out,"Compute[%3d]: %s, style = %s, group = %s\n",
i, compute[i]->id, compute[i]->style,
names[compute[i]->igroup]);
}
}
if (flags & DUMPS) {
int ndump = output->ndump;
Dump **dump = output->dump;
int *nevery = output->every_dump; \
char **vnames = output->var_dump;
char **names = group->names;
fprintf(out,"\nDump information:\n");
for (int i=0; i < ndump; ++i) {
fprintf(out,"Dump[%3d]: %s, file = %s, style = %s, group = %s, ",
i, dump[i]->id, dump[i]->filename,
dump[i]->style, names[dump[i]->igroup]);
if (nevery[i]) {
fprintf(out,"every = %d\n", nevery[i]);
} else {
fprintf(out,"every = %s\n", vnames[i]);
}
}
}
if (flags & FIXES) {
int nfix = modify->nfix;
Fix **fix = modify->fix;
char **names = group->names;
fprintf(out,"\nFix information:\n");
for (int i=0; i < nfix; ++i) {
fprintf(out,"Fix[%3d]: %s, style = %s, group = %s\n",
i, fix[i]->id, fix[i]->style, names[fix[i]->igroup]);
}
}
if (flags & VARIABLES) {
int nvar = input->variable->nvar;
int *style = input->variable->style;
char **names = input->variable->names;
char ***data = input->variable->data;
fprintf(out,"\nVariable information:\n");
for (int i=0; i < nvar; ++i) {
int ndata = 1;
fprintf(out,"Variable[%3d]: %-10s, style = %-10s, def =",
i,names[i],varstyles[style[i]]);
if ((style[i] != LOOP) && (style[i] != ULOOP))
ndata = input->variable->num[i];
for (int j=0; j < ndata; ++j)
fprintf(out," %s",data[i][j]);
fputs("\n",out);
}
}
if (flags & TIME) {
double wallclock = MPI_Wtime() - lmp->initclock;
double cpuclock = 0.0;
#if defined(_WIN32)
// from MSD docs.
FILETIME ct,et,kt,ut;
union { FILETIME ft; uint64_t ui; } cpu;
if (GetProcessTimes(GetCurrentProcess(),&ct,&et,&kt,&ut)) {
cpu.ft = ut;
cpuclock = cpu.ui * 0.0000001;
}
#else /* POSIX */
struct rusage ru;
if (getrusage(RUSAGE_SELF, &ru) == 0) {
cpuclock = (double) ru.ru_utime.tv_sec;
cpuclock += (double) ru.ru_utime.tv_usec * 0.000001;
}
#endif /* ! _WIN32 */
int cpuh,cpum,cpus,wallh,wallm,walls;
cpus = fmod(cpuclock,60.0);
cpuclock = (cpuclock - cpus) / 60.0;
cpum = fmod(cpuclock,60.0);
cpuh = (cpuclock - cpum) / 60.0;
walls = fmod(wallclock,60.0);
wallclock = (wallclock - walls) / 60.0;
wallm = fmod(wallclock,60.0);
wallh = (wallclock - wallm) / 60.0;
fprintf(out,"\nTotal time information (MPI rank 0):\n"
" CPU time: %4d:%02d:%02d\n"
" Wall time: %4d:%02d:%02d\n",
cpuh,cpum,cpus,wallh,wallm,walls);
}
if (flags & STYLES) {
available_styles(out, flags);
}
fputs("\nInfo-Info-Info-Info-Info-Info-Info-Info-Info-Info-Info\n\n",out);
// close output file pointer if opened locally thus forcing a hard sync.
if ((out != screen) && (out != logfile))
fclose(out);
}
void Info::available_styles(FILE * out, int flags)
{
fprintf(out,"\nStyles information:\n");
if(flags & ATOM_STYLES) atom_styles(out);
if(flags & INTEGRATE_STYLES) integrate_styles(out);
if(flags & MINIMIZE_STYLES) minimize_styles(out);
if(flags & PAIR_STYLES) pair_styles(out);
if(flags & BOND_STYLES) bond_styles(out);
if(flags & ANGLE_STYLES) angle_styles(out);
if(flags & DIHEDRAL_STYLES) dihedral_styles(out);
if(flags & IMPROPER_STYLES) improper_styles(out);
if(flags & KSPACE_STYLES) kspace_styles(out);
if(flags & FIX_STYLES) fix_styles(out);
if(flags & COMPUTE_STYLES) compute_styles(out);
if(flags & REGION_STYLES) region_styles(out);
if(flags & DUMP_STYLES) dump_styles(out);
if(flags & COMMAND_STYLES) command_styles(out);
}
void Info::atom_styles(FILE * out)
{
fprintf(out, "\nAtom styles:\n");
vector<string> styles;
for(Atom::AtomVecCreatorMap::iterator it = atom->avec_map->begin(); it != atom->avec_map->end(); ++it) {
styles.push_back(it->first);
}
print_columns(out, styles);
fprintf(out, "\n\n\n");
}
void Info::integrate_styles(FILE * out)
{
fprintf(out, "\nIntegrate styles:\n");
vector<string> styles;
for(Update::IntegrateCreatorMap::iterator it = update->integrate_map->begin(); it != update->integrate_map->end(); ++it) {
styles.push_back(it->first);
}
print_columns(out, styles);
fprintf(out, "\n\n\n");
}
void Info::minimize_styles(FILE * out)
{
fprintf(out, "\nMinimize styles:\n");
vector<string> styles;
for(Update::MinimizeCreatorMap::iterator it = update->minimize_map->begin(); it != update->minimize_map->end(); ++it) {
styles.push_back(it->first);
}
print_columns(out, styles);
fprintf(out, "\n\n\n");
}
void Info::pair_styles(FILE * out)
{
fprintf(out, "\nPair styles:\n");
vector<string> styles;
for(Force::PairCreatorMap::iterator it = force->pair_map->begin(); it != force->pair_map->end(); ++it) {
styles.push_back(it->first);
}
print_columns(out, styles);
fprintf(out, "\n\n\n");
}
void Info::bond_styles(FILE * out)
{
fprintf(out, "\nBond styles:\n");
vector<string> styles;
for(Force::BondCreatorMap::iterator it = force->bond_map->begin(); it != force->bond_map->end(); ++it) {
styles.push_back(it->first);
}
print_columns(out, styles);
fprintf(out, "\n\n\n");
}
void Info::angle_styles(FILE * out)
{
fprintf(out, "\nAngle styles:\n");
vector<string> styles;
for(Force::AngleCreatorMap::iterator it = force->angle_map->begin(); it != force->angle_map->end(); ++it) {
styles.push_back(it->first);
}
print_columns(out, styles);
fprintf(out, "\n\n\n");
}
void Info::dihedral_styles(FILE * out)
{
fprintf(out, "\nDihedral styles:\n");
vector<string> styles;
for(Force::DihedralCreatorMap::iterator it = force->dihedral_map->begin(); it != force->dihedral_map->end(); ++it) {
styles.push_back(it->first);
}
print_columns(out, styles);
fprintf(out, "\n\n\n");
}
void Info::improper_styles(FILE * out)
{
fprintf(out, "\nImproper styles:\n");
vector<string> styles;
for(Force::ImproperCreatorMap::iterator it = force->improper_map->begin(); it != force->improper_map->end(); ++it) {
styles.push_back(it->first);
}
print_columns(out, styles);
fprintf(out, "\n\n\n");
}
void Info::kspace_styles(FILE * out)
{
fprintf(out, "\nKSpace styles:\n");
vector<string> styles;
for(Force::KSpaceCreatorMap::iterator it = force->kspace_map->begin(); it != force->kspace_map->end(); ++it) {
styles.push_back(it->first);
}
print_columns(out, styles);
fprintf(out, "\n\n\n");
}
void Info::fix_styles(FILE * out)
{
fprintf(out, "\nFix styles:\n");
vector<string> styles;
for(Modify::FixCreatorMap::iterator it = modify->fix_map->begin(); it != modify->fix_map->end(); ++it) {
styles.push_back(it->first);
}
print_columns(out, styles);
fprintf(out, "\n\n\n");
}
void Info::compute_styles(FILE * out)
{
fprintf(out, "\nCompute styles:\n");
vector<string> styles;
for(Modify::ComputeCreatorMap::iterator it = modify->compute_map->begin(); it != modify->compute_map->end(); ++it) {
styles.push_back(it->first);
}
print_columns(out, styles);
fprintf(out, "\n\n\n");
}
void Info::region_styles(FILE * out)
{
fprintf(out, "\nRegion styles:\n");
vector<string> styles;
-#define REGION_CLASS
-#define RegionStyle(key,Class) styles.push_back(#key);
-#include "style_region.h"
-#undef REGION_CLASS
+
+ for(Domain::RegionCreatorMap::iterator it = domain->region_map->begin(); it != domain->region_map->end(); ++it) {
+ styles.push_back(it->first);
+ }
+
print_columns(out, styles);
fprintf(out, "\n\n\n");
}
void Info::dump_styles(FILE * out)
{
fprintf(out, "\nDump styles:\n");
vector<string> styles;
for(Output::DumpCreatorMap::iterator it = output->dump_map->begin(); it != output->dump_map->end(); ++it) {
styles.push_back(it->first);
}
print_columns(out, styles);
fprintf(out, "\n\n\n");
}
void Info::command_styles(FILE * out)
{
fprintf(out, "\nCommand styles (add-on input script commands):\n");
vector<string> styles;
for(Input::CommandCreatorMap::iterator it = input->command_map->begin(); it != input->command_map->end(); ++it) {
styles.push_back(it->first);
}
print_columns(out, styles);
fprintf(out, "\n\n\n");
}
/* ---------------------------------------------------------------------- */
// the is_active() function returns true if the selected style or name
// in the selected category is currently in use.
bool Info::is_active(const char *category, const char *name)
{
if ((category == NULL) || (name == NULL)) return false;
const char *style = "none";
const int len = strlen(name);
if (strcmp(category,"package") == 0) {
if (strcmp(name,"gpu") == 0) {
return (modify->find_fix("package_gpu") >= 0) ? true : false;
} else if (strcmp(name,"intel") == 0) {
return (modify->find_fix("package_intel") >= 0) ? true : false;
} else if (strcmp(name,"kokkos") == 0) {
return (lmp->kokkos && lmp->kokkos->kokkos_exists) ? true : false;
} else if (strcmp(name,"omp") == 0) {
return (modify->find_fix("package_omp") >= 0) ? true : false;
} else error->all(FLERR,"Unknown name for info package category");
} else if (strcmp(category,"newton") == 0) {
if (strcmp(name,"pair") == 0) return (force->newton_pair != 0);
else if (strcmp(name,"bond") == 0) return (force->newton_bond != 0);
else if (strcmp(name,"any") == 0) return (force->newton != 0);
else error->all(FLERR,"Unknown name for info newton category");
} else if (strcmp(category,"pair") == 0) {
if (force->pair == NULL) return false;
if (strcmp(name,"single") == 0) return (force->pair->single_enable != 0);
else if (strcmp(name,"respa") == 0) return (force->pair->respa_enable != 0);
else if (strcmp(name,"manybody") == 0) return (force->pair->manybody_flag != 0);
else if (strcmp(name,"tail") == 0) return (force->pair->tail_flag != 0);
else if (strcmp(name,"shift") == 0) return (force->pair->offset_flag != 0);
else error->all(FLERR,"Unknown name for info pair category");
} else if (strcmp(category,"comm_style") == 0) {
style = commstyles[comm->style];
} else if (strcmp(category,"min_style") == 0) {
style = update->minimize_style;
} else if (strcmp(category,"run_style") == 0) {
style = update->integrate_style;
} else if (strcmp(category,"atom_style") == 0) {
style = atom->atom_style;
} else if (strcmp(category,"pair_style") == 0) {
style = force->pair_style;
} else if (strcmp(category,"bond_style") == 0) {
style = force->bond_style;
} else if (strcmp(category,"angle_style") == 0) {
style = force->angle_style;
} else if (strcmp(category,"dihedral_style") == 0) {
style = force->dihedral_style;
} else if (strcmp(category,"improper_style") == 0) {
style = force->improper_style;
} else if (strcmp(category,"kspace_style") == 0) {
style = force->kspace_style;
} else error->all(FLERR,"Unknown category for info is_active()");
int match = 0;
if (strcmp(style,name) == 0) match = 1;
if (!match && lmp->suffix_enable) {
if (lmp->suffix) {
char *name_w_suffix = new char [len + 2 + strlen(lmp->suffix)];
sprintf(name_w_suffix,"%s/%s",name,lmp->suffix);
if (strcmp(style,name_w_suffix) == 0) match = 1;
delete[] name_w_suffix;
}
if (!match && lmp->suffix2) {
char *name_w_suffix = new char [len + 2 + strlen(lmp->suffix2)];
sprintf(name_w_suffix,"%s/%s",name,lmp->suffix2);
if (strcmp(style,name_w_suffix) == 0) match = 1;
delete[] name_w_suffix;
}
}
return match ? true : false;
}
/* ---------------------------------------------------------------------- */
// the is_available() function returns true if the selected style
// or name in the selected category is available for use (but need
// not be currently active).
bool Info::is_available(const char *category, const char *name)
{
if ((category == NULL) || (name == NULL)) return false;
const int len = strlen(name);
int match = 0;
if (strcmp(category,"command") == 0) {
if (input->command_map->find(name) != input->command_map->end())
match = 1;
} else if (strcmp(category,"compute") == 0) {
if (modify->compute_map->find(name) != modify->compute_map->end())
match = 1;
if (!match && lmp->suffix_enable) {
if (lmp->suffix) {
char *name_w_suffix = new char [len + 2 + strlen(lmp->suffix)];
sprintf(name_w_suffix,"%s/%s",name,lmp->suffix);
if (modify->compute_map->find(name_w_suffix) != modify->compute_map->end())
match = 1;
delete[] name_w_suffix;
}
if (!match && lmp->suffix2) {
char *name_w_suffix = new char [len + 2 + strlen(lmp->suffix2)];
sprintf(name_w_suffix,"%s/%s",name,lmp->suffix2);
if (modify->compute_map->find(name_w_suffix) != modify->compute_map->end())
match = 1;
delete[] name_w_suffix;
}
}
} else if (strcmp(category,"fix") == 0) {
if (modify->fix_map->find(name) != modify->fix_map->end())
match = 1;
if (!match && lmp->suffix_enable) {
if (lmp->suffix) {
char *name_w_suffix = new char [len + 2 + strlen(lmp->suffix)];
sprintf(name_w_suffix,"%s/%s",name,lmp->suffix);
if (modify->fix_map->find(name_w_suffix) != modify->fix_map->end())
match = 1;
delete[] name_w_suffix;
}
if (!match && lmp->suffix2) {
char *name_w_suffix = new char [len + 2 + strlen(lmp->suffix2)];
sprintf(name_w_suffix,"%s/%s",name,lmp->suffix2);
if (modify->fix_map->find(name_w_suffix) != modify->fix_map->end())
match = 1;
delete[] name_w_suffix;
}
}
} else if (strcmp(category,"pair_style") == 0) {
if (force->pair_map->find(name) != force->pair_map->end())
match = 1;
if (!match && lmp->suffix_enable) {
if (lmp->suffix) {
char *name_w_suffix = new char [len + 2 + strlen(lmp->suffix)];
sprintf(name_w_suffix,"%s/%s",name,lmp->suffix);
if (force->pair_map->find(name_w_suffix) != force->pair_map->end())
match = 1;
delete[] name_w_suffix;
}
if (!match && lmp->suffix2) {
char *name_w_suffix = new char [len + 2 + strlen(lmp->suffix2)];
sprintf(name_w_suffix,"%s/%s",name,lmp->suffix2);
if (force->pair_map->find(name_w_suffix) != force->pair_map->end())
match = 1;
delete[] name_w_suffix;
}
}
} else error->all(FLERR,"Unknown category for info is_available()");
return match ? true : false;
}
/* ---------------------------------------------------------------------- */
// the is_defined() function returns true if a particular ID of the
// selected category (e.g. fix ID, group ID, region ID etc.) has been
// defined and thus can be accessed. It does *NOT* check whether a
// particular ID has a particular style.
bool Info::is_defined(const char *category, const char *name)
{
if ((category == NULL) || (name == NULL)) return false;
if (strcmp(category,"compute") == 0) {
int ncompute = modify->ncompute;
Compute **compute = modify->compute;
for (int i=0; i < ncompute; ++i) {
if (strcmp(compute[i]->id,name) == 0)
return true;
}
} else if (strcmp(category,"dump") == 0) {
int ndump = output->ndump;
Dump **dump = output->dump;
for (int i=0; i < ndump; ++i) {
if (strcmp(dump[i]->id,name) == 0)
return true;
}
} else if (strcmp(category,"fix") == 0) {
int nfix = modify->nfix;
Fix **fix = modify->fix;
for (int i=0; i < nfix; ++i) {
if (strcmp(fix[i]->id,name) == 0)
return true;
}
} else if (strcmp(category,"group") == 0) {
int ngroup = group->ngroup;
char **names = group->names;
for (int i=0; i < ngroup; ++i) {
if (strcmp(names[i],name) == 0)
return true;
}
} else if (strcmp(category,"region") == 0) {
int nreg = domain->nregion;
Region **regs = domain->regions;
for (int i=0; i < nreg; ++i) {
if (strcmp(regs[i]->id,name) == 0)
return true;
}
} else if (strcmp(category,"variable") == 0) {
int nvar = input->variable->nvar;
char **names = input->variable->names;
for (int i=0; i < nvar; ++i) {
if (strcmp(names[i],name) == 0)
return true;
}
} else error->all(FLERR,"Unknown category for info is_defined()");
return false;
}
static void print_columns(FILE* fp, vector<string> & styles)
{
if (styles.size() == 0) {
fprintf(fp, "\nNone");
return;
}
std::sort(styles.begin(), styles.end());
int pos = 80;
for (int i = 0; i < styles.size(); ++i) {
// skip "secret" styles
if (isupper(styles[i][0])) continue;
int len = styles[i].length();
if (pos + len > 80) {
fprintf(fp,"\n");
pos = 0;
}
if (len < 16) {
fprintf(fp,"%-16s",styles[i].c_str());
pos += 16;
} else if (len < 32) {
fprintf(fp,"%-32s",styles[i].c_str());
pos += 32;
} else if (len < 48) {
fprintf(fp,"%-48s",styles[i].c_str());
pos += 48;
} else if (len < 64) {
fprintf(fp,"%-64s",styles[i].c_str());
pos += 64;
} else {
fprintf(fp,"%-80s",styles[i].c_str());
pos += 80;
}
}
}